location measure로 outlier에 robust한 median(or weighted median)이나 trimmed mean을 mean대신 사용해야 할 수도 있다.
mean(state[["Population"]])## [1] 6162876
mean(state[["Population"]], trim = 0.1)## [1] 4783697
median(state[["Population"]])## [1] 4436370
# install.packages("matrixStats")
weighted.mean(state[["Murder.Rate"]], w = state[["Population"]])## [1] 4.445834
matrixStats::weightedMedian(state[["Murder.Rate"]], w = state[["Population"]])## [1] 4.4
분산과 표준편차는 변동성 통계량중에서 가장 넓게 퍼져있고 outlier에 민감하다.
sd(state[["Population"]])## [1] 6848235
IQR(state[["Population"]])## [1] 4847308
mad(state[["Population"]])## [1] 3849870
quantile(state[["Murder.Rate"]], probs = c(.05, .25, .5, .75, .95))## 5% 25% 50% 75% 95%
## 1.600 2.425 4.000 5.550 6.510
boxplot(state[["Population"]] / 1000000, ylab = "Population (millions)")breaks <- seq(min(state[["Population"]]), max(state[["Population"]]), length.out = 11)
pop_freq <- cut(state[["Population"]], breaks = breaks, right = TRUE, include.lowest = TRUE)
table(pop_freq) # bin의 크기가 첫번째 bin을 제외하고 37로 동일## pop_freq
## [5.64e+05,4.23e+06] (4.23e+06,7.9e+06] (7.9e+06,1.16e+07]
## 24 14 6
## (1.16e+07,1.52e+07] (1.52e+07,1.89e+07] (1.89e+07,2.26e+07]
## 2 1 1
## (2.26e+07,2.62e+07] (2.62e+07,2.99e+07] (2.99e+07,3.36e+07]
## 1 0 0
## (3.36e+07,3.73e+07]
## 1
빈도표와 백분위수 모두 데이터를 bin을 생성함으로써 요약한다. 차이점은 빈도표는 bin의 크기가 모두 동일한 반면 백분위수는 bin의 크기가 다르고 bin마다 빈도 수가 동일하다.
hist(state[["Population"]], breaks = breaks)hist(state[["Murder.Rate"]], freq = FALSE)
lines(density(state[["Murder.Rate"]]), lew = 3, col = "blue")## Warning in plot.xy(xy.coords(x, y), type = type, ...): "lew"는 그래픽 매개변
## 수가 아닙니다
Converting numeric data to categorical data is an important and widely used step in data analysis since it reduces the complexity (and size) of the data
dfw_airline %>%
gather(cause, val) %>%
ggplot(aes(cause, val)) +
geom_col() +
scale_y_continuous(
labels = function(y) y / 6 # 6년간 데이터여서
)etfs <- sp500_px[row.names(sp500_px)>"2012-07-01",
sp500_sym[sp500_sym$sector=="etf", 'symbol']]
library(corrplot)
corrplot(cor(etfs), method = "ellipse")cor함수도 trim인자를 제공한다.
데이터가 작다면 이상치에 robust한 Spearman’s rho or Kendall’s tau도 있다.
주의: 데이터를 랜덤하게 재배치하면 우연히 상관계수가 양수나 음수가 될 수도 있다.
아주 큰 데이터에서는 산점도 보다는..
kc_tax0 <- kc_tax %>%
filter(TaxAssessedValue < 750000 & SqFtTotLiving>100 & SqFtTotLiving<3500)
ggplot(kc_tax0, aes(SqFtTotLiving, TaxAssessedValue)) +
geom_hex(color = "white") +
theme_bw() +
# scale_fill_viridis_c()
scale_fill_gradient(low = "white", high = "blue") +
labs(x = "Finished Squared Feet", y = "Tax Assessed Value")## Warning: package 'hexbin' was built under R version 3.5.2
main 구름위에 두번째 구름이 집의 크기는 동일하지만 더 높은 과세 평가액을 갖는다.
ggplot(kc_tax0, aes(SqFtTotLiving, TaxAssessedValue)) +
theme_bw() +
geom_point(alpha = 0.1) +
geom_density2d(color = "white") +
labs(x = "Finished Square Feet", y = "Tax Assessed Value") 두번째 peak를 주목하자!
library(descr)## Warning: package 'descr' was built under R version 3.5.2
x_tab <- CrossTable(lc_loans$grade, lc_loans$status, prop.c = FALSE, prop.chisq = FALSE, prop.t = FALSE)ggplot(airline_stats, aes(airline, pct_carrier_delay)) +
geom_boxplot()## Warning: Removed 28 rows containing non-finite values (stat_boxplot).
kc_tax0 %>%
filter(ZipCode %in% c(98188, 98105, 98108, 98126)) %>%
ggplot(aes(SqFtTotLiving, TaxAssessedValue)) +
geom_hex(color = "white") +
theme_bw() +
scale_fill_gradient(low = "white", high = "blue") +
labs(x="Finished Square Feet", y="Tax Assessed Value") +
facet_wrap("ZipCode")sample bias를 조심하라.
random chance에 의한 error와 bias에 의한 error의 구분은 아주 중요하다.
유동적인 상황에서 표본추출을 해야 할 경우 시기가 중요할 수 있다.(평일 오전 10시의 웹 방문자와 주말 오후 10시의 웹 방문자가 다를 수 있다.)
빅데이터 시대라고 해도 의외로 데이터 개수가 적을수록 더 유리한 경우가 있다. 랜덤표본추출에 시간과 노력을 기울일수록 편향이 줄 뿐만 아니라 데이터 탐색 및 게이터 품질에 더 집중할 수 있다.
빅데이터의 가치는 데이터가 크고 동시에 희박할 때이다.
selection bias의 종류 - vast search effect: 큰 데이터 셋에서 반복적으로 다른 모델을 만들고 다른 질문을 하면 흥미로운 것을 찾기 마련이다. 의미가 있는 것일까? 우연히 얻은 것일까? - nonrandom sampling(self-selection sampling bias) - cherry-picking data - selection of time intervals that accentuate a particular statistical effect - stopping an experiment when the results look “interesting”
해결 방안: holdout set, permutation test
평균으로의 회귀는 어떤 변수를 연속적으로 측정했을 때 나타나는 현상이다. 예외 경우를 특별히 생각하고 의미를 부여하는 것은 선택편향으로 이어진다.
평균으로의 회귀는 선택편향의 특정한 형태의 결과이다.
우리의 추정치나 모델은 표본에 기반하기 때문에 에러가 있을 수 있다. 우리가 다른 표본을 뽑았다면 결과가 다를 수 있다. 그래서 주된 관심사는 표본 변동성(sampling variability)이다.
# samp_mean_05 <- loans_income %>%
# sample_n(5000) %>%
# mutate(idx = rep(1:1000, each = 5)) %>%
# group_by(idx) %>%
# summarise(mean(x))
loans_income <- loans_income$x
# take a simple random sample
samp_data <- vector("list", 3)
samp_data[[1]] <- tibble(
income = sample(loans_income, 1000),
type = "data_dist"
)
samp_data[[2]] <- tibble(
income = tapply(sample(loans_income, 1000 * 5), rep(1:1000, each = 5), mean),
type = "mean_of_5")
samp_data[[3]] <- tibble(
income = tapply(sample(loans_income, 1000 * 20), rep(1:1000, each = 20), mean),
type = "mean_of_20")
income <- bind_rows(samp_data) %>%
mutate(type = factor(type,
levels = c("data_dist", "mean_of_5", "mean_of_20"),
labels = c("Data", "Mean of 5", "Mean of 20")))
ggplot(income, aes(income)) +
geom_histogram(bins = 40) +
facet_wrap(~type, ncol = 1) +
theme_bw()가설검정과 신뢰구간에 있어 CLT의 중요성을 인지하고 있어야 한다. 하지만 표준오차를 추정함에 있어 bootstrap을 사용하면 된다.
library(boot)
stat_fun <- function(x, idx) median(x[idx])
boot_obj <- boot(loans_income, R = 1000, statistic = stat_fun)bootstrap의 주의할 점은 최초로 주어진 데이터가 모분포에 근사해야한다는 것이다. 참고로 bias가 양수면 과대추정한다.
모델 파라미터의 안정성(변동성)을 추정하거나 예측력을 높이기 위해, bootstrap데이터를 가지고 모델을 돌려볼 수 있다. 분류 및 회귀 트리를 사용할 때, 여러 부트스트랩 샘플을 가지고 트리를 여러 개 만든 다음 각 트리에서 나온 예측값을 평균 내는 것이(분류 문제에서는 과반수 투표를 한다) 일반적으로 단일 트리를 사용하는 것보다 효과적이다. 이 프로세스를 bagging이라고 한다.
주의: 부트스트랩은 표본크기가 작은 것을 보완하기 위한 것이 아니다. 새 데이터를 만드는 것도 아니며 기존 데이터 집합의 빈 곳을 채우는 것도 아니다. 모집단에서 추가적으로 표본을 뽑는다고 할 때, 그 표본이 얼마나 원래 표본과 비슷할지를 알려줄 뿐이다.
부트스트랩 신뢰구간 구하기 1. Draw a radom sample of size n with replacement from the data 2. Record the statistic of interest for the resample 3. Repeat steps 1-2 many(R) times. 4. For an x% CI, trim[(1 - [x/100]) / 2]% of the R resample results from either end of the distribution. 5. The trim points are the endpoints of an x% bootstrap CI
분포가 소득 데이터와 같이 비대칭하게 기울어져 있거나 이항 데이터같이 이산적일 수 있다. 대칭 및 비대칭 분포 모두 긴 꼬리를 가질 수 있다. Nassim Taleb는 주식시장의 붕괴와 같은 이례적인 사건이, 정규분포로 예측되는 것보다 훨씬 더 자주 일어날 수 있다고 예측하는 black swan theory를 제안했다.
주가 수익률은 데이터의 긴 꼬리 특성을 설명하기 위한 종은 예가 된다.
nflx <- sp500_data[, "NFLX"]$NFLX
nflx <- diff(log(nflx[nflx > 0]))
qqnorm(nflx)
abline(0, 1, col = "grey")
Figure 1: QQ-Plot of the returns for NFLX
관측된 자료에 적합한 통계 분포를 찾는 작업에 관한 통계적 문헌이 많이 있으므로 너무 데이터만을 보고 판단하지 말자. 예를 들어, 5초마다 연속적으로 서버의 인터넷 트래픽 수준에 대한 데이터를 얻는다고 하자. 포아송 분포가 ’events per time period’를 모델링하는 데 가장 적합한 분포라는 사전 지식이 있다면 큰 도움이 될 것이다.
사건 발생률은 시간에 따라 일정하지 않다. 하지만 변화 주기가 일반적인 사건 발생 구간보다 훨씬 길다면 문제가 안 된다. 앞서 언급했듯이 비율이 상대적으로 일정한 구간으로 분석을 세분화하면 되기 때문이다. 그러나 사건 발생률이 시간에 따라 지속적으로 변한다면 지수 (또는 포아송)분포는 유용하지 않다. 기계 고장이 대표적이다. 시간이 지날수록 고장 위험은 증가한다. weibull dist는 지수분포를 확장한 것으로, shape parameter \(\beta\)로 발생률이 달라질 수 있다. \(\beta\) > 1 일 경우, 발생률은 시간이 지남에 따라 증가하며, \(\beta\) < 1이면 감소한다. weibull 분포는 사건 발생률 대신 고장 시간 분석에 사용되기 때문에 두 번째 인자는 구간당 사건 발생률보다는 특성 수명으로 표현된다.
A/B Test는 두 개의 treatments, 제품, 절차 중 어느 쪽이 더 우월한 것을 입증하기 위해 두 개의 그룹으로 나누어 진행하는 실험이다.
Revenue/page-view with price A: mean = 3.87, SD = 51.10 Revenue/page-view with price B: mean = 4.11, SD = 62.98
값이 음수가 될 수 없는데 음수가 될 수 있다고 제안하고 있으므로 표준편차는 유용하지 않다. 이 경우에는 mean absolute deviation이 더 합리적이다.
단일 지표 또는 검정통계량을 사전에 미리 정해놓아야 한다. 실험을 수행한 뒤 나중에 검정통계량을 선택한다면 연구자 편향이라는 함정에 빠지게 된다.
A/B 검정은 마케팅 및 전자 상거래 분야에서 널리 사용되고 있지만 그렇다고 이것이 유일한 통계 실험 유형인 것은 아니다. 추가적인 treatment를 포함할 수 있고, 반복 측정을 할 수 있다. 제약회사의 임상 실험과 같이 대상이 매우 귀하고 비용이 비싸며 측정에 많은 시간이 필요할 경우, 실험을 중간에 중단하고 결론을 얻을 수 있는 장치를 마련해두고 실험을 설계한다.
데이터 과학자는 ‘가격 A와 가격 B의 차이가 통계적으로 유의한가?’ 보다는 ’가능한 여러 가격 중에서 가장 좋은 가격은 얼마일까?’에 관심이 있다. 이러한 실험을 위해서는 Multi-Arm Bandit Algorithm이 필요
사람은 랜덤성을 과소평가하려고 한다. 예를 들어, 동전 던지기에서 앞면이 6번 연속으로 나오면 뭔가 의미 있고 단순한 우연은 아닐 것이다라고 생각하는 경향이 있다는 점이다.
통계학에서 resampling이란 랜덤한 변동성을 알아보자는 일반적인 목표를 가지고, 관찰된 데이터의 값에서 표본을 반복적으로 추출하는 것을 의미한다. 또한 일부 머신러닝 모델의 정확성을 평가하고 향상시키는 데에도 적용할 수 있다. resamling에는 크게 bootstrap과 permutation test가 있다. bootstrap은 추정치의 신뢰도를 평가하기 위해 이요된다. permutation test는 일반적으로 두 개 이상의 그룹과 관련된 가설을 검증하는 데 사용된다.
permutation 과정에서는 두 개 이상의 표본이 관여되며 이들은 통상적으로 A/B 또는 기타 가설검정을 위해 사용되는 그룹들이다. 1. 여러 그룹의 결과를 단일 데이터 집합으로 결합한다. 2. 결합된 데이터를 잘 섞은 후, 그룹 A와 동일한 크기의 표본을 비복원 추출한다. 3. 나머지 데이터에서 그룹 B와 동일한 크기의 샘플을 비복원 추출한다. 4. C, D 등의 그룹에 대해서도 동일한 작업을 수행한다. 5. 원래 샘플에 대해 구한 통계량 또는 추정치가 무엇이었든 간에 지금 추출한 재표본에 대해 모두 다시 계산하고 기록한다. 이것으로 한 번의 permutation 반복이 진행된다. 6. 앞서 단계들을 R번 반복하여 검정통계량의 순열 분포를 얻는다.
관찰된 차이가 순열 분포 바깥에 있다면, 우연 때문이 아니라고 결론 내릴 수 있다: statistically significant
Example: Web Stickiness
landing page: 온라인광고에 노출된 사용자가 광고메시지에 관심을 가지고 보다 많은 정보 탐색 또는 광고주가 의도한 행동을 하기 위해 광고물을 클릭했을 때 처음 접하게 되는 화면을 랜딩페이지라고 한다.
랜딩페이지의 중요성: 1. 사용자의 의도에 맞는 컨텐츠를 강조하여 방문자에게 확실히 전달 가능 2. 관심 방문자의 상담 또는 문의 내용 DB 수집 용이 3. 유저들의 원하는 정보 및 관심사를 몇 번의 클릭 없이 바로 접근 가능 4. 상세 의도를 가지고 방문한 유저들의 이탈 최소화 기대
web_page_data <- web_page_data %>%
mutate(Time = 60 * Time)
ggplot(web_page_data, aes(Page, Time)) +
geom_boxplot()mean_a <- mean(unlist(web_page_data[web_page_data$Page == "Page A", "Time"]))
mean_b <- mean(unlist(web_page_data[web_page_data$Page == "Page B", "Time"]))
mean_b - mean_a## [1] 21.4
# permutation test
perm_fun <- function(x, n1, n2) {
n <- n1 + n2
idx_b <- sample(1:n, n1)
idx_a <- setdiff(1:n, idx_b)
mean_diff <- mean(x[idx_b]) - mean(x[idx_a])
return(mean_diff)
}
perm_diffs <- rep(0, 1000)
for(i in 1:1000)
perm_diffs[i] <- perm_fun(web_page_data$Time, 21, 15)
hist(perm_diffs, xlab='Session time differences (in seconds)')
abline(v = mean_b - mean_a)permutation test의 변형
exhaustive permutation test = exact test
통계적 유의성이란, 통계학자가 자신의 실험 결과가 우연히 일어난 것인지 아니면 우연히 일어날 수 없는 극단적인 건인지를 판단하는 방법이다. 결과가 우연히 벌어질 수 있는 변동성의 바깥에 존재한다면 우리는 이것을 통계적으로 유의하다고 말한다.
가격 A가 가격 B에 비해 0.0368%p 개선되었다. 전환율이 너무 낮아, 실제 필요한 표본크기를 결정하는 데 매우 중요한 값은 200개 정도에 불과하다. 두 전환율 사이에 차이가 있는지를 permutation test를 통해 확인해보자.
obs_pct_diff <- 100*(200/23739 - 182/22588)
conversion <- c(rep(0, 45945), rep(1, 382))
perm_diffs <- rep(0, 1000)
for(i in 1:1000)
perm_diffs[i] <- 100*perm_fun(conversion, 23739, 22588 )
hist(perm_diffs, xlab='Session time differences (in seconds)')
abline(v = obs_pct_diff)# p-value
mean(perm_diffs > obs_pct_diff)## [1] 0.307
하지만, 이 경우에는 이항분포를 따르기 때문에 정규근사를 이용할 수 있다.
prop.test(x = c(200, 182), n = c(23739, 22588), alternative = "greater")##
## 2-sample test for equality of proportions with continuity
## correction
##
## data: c(200, 182) out of c(23739, 22588)
## X-squared = 0.14893, df = 1, p-value = 0.3498
## alternative hypothesis: greater
## 95 percent confidence interval:
## -0.001057439 1.000000000
## sample estimates:
## prop 1 prop 2
## 0.008424955 0.008057376
p 값의 의미: 랜덤 모델이 주어졌을 때, 그 결과가 관찰된 결과보다 더 극단적일 확률
다중검정시 \(1 - (1 - \alpha)^n\) 문제에 주의하자. 이것은 오버피팅과 관련이 있다.
중복 문제:
false discovery rate은 원래 주어진 여러 개의 가설검정들 가운데 하나가 유의미한 효과가 있다고 잘못 판단하는 비율을 나타내는 데 사용되었다. 이것은 분류 문제를 다루는 데에도 사용된다. 이때의 false discovery은 한 레코드를 잘못 분류한 것을 의미한다.
library(lmPerm)## Warning: package 'lmPerm' was built under R version 3.5.2
summary(aovp(Time ~ Page, data = four_sessions))## [1] "Settings: unique SS "
## Component 1 :
## Df R Sum Sq R Mean Sq Iter Pr(Prob)
## Page1 3 831.4 277.13 5000 0.0804 .
## Residuals 16 1618.4 101.15
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
여기서 p 값은 재표집된 분산이 관찰된 데이터의 분산을 초과하는 비율이다.
summary(aov(Time ~ Page, data = four_sessions))## Df Sum Sq Mean Sq F value Pr(>F)
## Page 3 831.4 277.1 2.74 0.0776 .
## Residuals 16 1618.4 101.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ANOVA의 결과 중 유용한 점 중 하나는 그룹 처리, 상호작용 효과, 오차와 관련된 분산의 구성 요소들을 구분하는 데 있다.
clicks <- xtabs(Rate ~ Click + Headline, data = click_rates)
chisq.test(clicks, simulate.p.value = TRUE)##
## Pearson's Chi-squared test with simulated p-value (based on 2000
## replicates)
##
## data: clicks
## X-squared = 1.6659, df = NA, p-value = 0.4793
chisq.test(clicks)##
## Pearson's Chi-squared test
##
## data: clicks
## X-squared = 1.6659, df = 2, p-value = 0.4348
fisher.test(clicks)##
## Fisher's Exact Test for Count Data
##
## data: clicks
## p-value = 0.4824
## alternative hypothesis: two.sided
대부분의 실험에서 통계적 유의성을 조사하는 것이 목적이 아니라 최적의 treatment를 찾는 것이다. multi-armed bandits이 더 정확한 해결책이다. 카이제곱 검정의 데이터사이언스로의 응용은 특히 피셔의 정확검정을 활용하는 예로 웹 실험에 적합한 표본크기를 판별하는 일을 들 수 있다. 이러한 실험은 종종 클릭율이 매우 낮기 때문에 수천 번의 싫머에도 불구하고 집계 비율이 너무 낮아 실험을 통해 확실한 결론을 내리기 어렵다. 이러한 경우 피셔의 정확검정, 카이제곱검정, 그리고 기타 검정은 검정력이나 표본크기를 계산하는 데 유용할 수 있다.
MAB알고리즘은 실험설계에 대한 전통적인 통계적 접근 방식보다 명시적인 최적화와 좀 더 빠른 의사 결정을 가능하게 하며, 여러 테스트, 특히 웹 테스트를 위해 이를 사용한다. epsilon-greedy algorithm과 좀 더 복잡한 알고리즘인 Thompson’s samling을 사용하는 알고리즘이 있다. 톰슨 샘플링은 베이지안 방식을 사용한다. beta 분포를 사용하여 수익의 일부 사전 분포를 가정한다. 각 추출 정보가 누적되면서 정보가 업데이트되기 때문에, 다음번에 최고 손잡이를 선택할 확률을 효과적으로 최적화할 수 있다.
검정력이란 바로 특정 표본 조건에서 특저안 효과크기를 알아낼 수 있는 확률을 의미한다. A/B test을 위해 데이터를 수집하고 처리하는 데 비용이 발생하는 경우 사용한다. 1. 최대한 (사전 정보를 이용해서) 결과 데이터가 비슷하게 나올 수 있는 가상의 데이터를 생각해보자. 예를 들면 타자를 위해 20개의 1과 80개의 0이 들어 있는 상자를 생각한다든지, 아니면 웹 페이제 방문 시간을 관측한 자료가 담겨 있는 상자를 생각할 수 있다. 2. 첫 표본에서 원하는 효과크기를 더해서 두 번째 표본을 만든다. 예를 들면 33개의 1과 67개의 0을 가진 두 번째 상자, 혹은 각 초기 방문시간에 25초를 더한 두 번째 상자를 만들 수 있다. 3. 각 상자에서 크기 n인 부트스트렙 표본을 추출한다 4. 두 부트스트렙 표본에 대해서 순열 가설검정을 진행한다(혹은 수식 기반의 가설검정). 그리고 여기에 통계적으로 유의미한 차이가 있는지 기록한다. 5. 3~4단계를 여러 번 반복한 후, 얼마나 자주 유의미한 차이가 발견되는지 알아본다. 이 확률이 바로 검정력 추정치다.
# 표본크기 정하기
pwr.2p.test(h = ..., n = ..., sig.level = ..., power = )
h= effect size (as a proportion)
n = sample size
sig.level = the significance level (alpha) at which the test will be conducted
power = power (probability of detecting the effect size)Cross-Validation
penalized regression는 개념적으로 AIC와 같다. 개별 모델 집합들을 명시적으로 검색하는 대신 모델 적합 방정식에 많은 변수에 대해 모델에 불이익을 주는 제약 조건을 추가한다. 단계적, 전진선택, 후진선택처럼 예측변수를 완전히 제거하는 대신, 벌정회귀에서는 계수 크기를 감소시키거나 경우에 따라 거의 0으로 만들어 벌점을 적용한다. ridge regression and lasso regression.
stepwise and all subset regression are in-sample methods to assess and tune models. 오버피팅 될 수 있으며, 새 데이터를 적용할 때 잘 맞지 않을 수도 있다. In linear regression, overfitting is typically not a major issue, due to the simple (linear) global structure imposed on the data. For more sophisticated types of models, particularly iterative procedures that respond to local data structure, CV is a very important tool.
Weighted Regression
복잡한 설문 분석에 중요하다.
library(lubridate)
house$year <- year(house$DocumentDate)
# 2005년 이래 지난 연수를 가중치로 사용
house$Weight <- house$Year - 2005
house_wt <- lm(AdjSalePrice ~ SqFtTotLiving + SqFtLot + Bathrooms + Bedrooms + BldgGrade, data=house, weight=Weight)
round(cbind(house_lm=house_lm$coefficients,
house_wt=house_wt$coefficients), digits=3)Reference coding: 통계학자들이 많이 사용하는 코딩 형태. 여기서 한 요인을 기준으로 하고 다른 요인들이 이 기준에 따라 비교할 수 있도록 한다(유의어: treatment coding). 첫 번째 팩터 레벨을 reference로 사용한다. One hot encoder: 머신러닝 분야에서 많이 사용되는 코딩. 모든 요인 수준이 계속 유지된다. 특정 머신러닝 알고리즘에는 유용한 반면, 다중회귀분석에는 적절하지 않다. nearest neighbor와 tree model에서 사용 Deviation coding: 기준 수준과는 반대로 전체 평균에 대해 각 수준을 비교하는 부호화 방법 Polynomial coding: ordered factor에 적절하다.
Factor Variables with Many Levels
어떤 요인변수는 가능한 수준의 수가 많아, 많은 양의 이진 더미를 생성할 수 있다. 예를 들면, 우편번호는 요인변수이며 미국에는 43000개의 우편번호가 있다. 이러한 경우 데이터와 예측변수와 결과 간의 관계를 탐색하여 유용한 정보가 범주에 포함되는지 여부를 판단하는 것이 유용하다. 그렇다면 모든 요소를 유지하는 것이 좋을지 아니면 수준을 통합하는 것이 나을지를 결정해야 한다.
ZipCode는 주택 가격에 대한 위치의 효과의 proxy이기 때문에 중요한 변수이다. 우편번호 중 몇몇은 매물이 하나뿐인 경우도 있다. 따라서 통합이 필요하다. 대도시 지역의 지리적 위치에 대응되는 처음 두 자리 또는 세 자리만을 사용하여 우편번호를 통합할 수 있다. 하지만 King County에 경우 거의 모든 판매가 980xx 또는 981xx에서 이루어지기 때문에 이런 방법이 큰 도움이 되지 않는다.
대안은 매매 가격과 같은 다른 변수에 따라 우편번호를 그룹으로 묶는 것이다. 아니면 초기 모델의 잔차를 사용하여 우편번호 그룹을 만드는 방법도 좋다.
zip_groups <- house %>%
mutate(resid = residuals(house_lm)) %>%
group_by(ZipCode) %>%
summarise(med_resid = median(resid), # 이상치에 덜 민감한 중앙값
cnt = n()) %>%
arrange(med_resid) %>%
mutate(cum_cnt = cumsum(cnt),
ZipGroup = ntile(cum_cnt, 5))
house <- house %>%
left_join(select(zip_groups, ZipCode, ZipGroup), by = "ZipCode")Bedrooms의 계수가 음수이다. 침실 계수를 늘릴수록 그 가치가 감소한다는 것을 의미한다. 이는 예측변수들이 서로 연관되어 있기 때문이다. 집이 클수록 침실이 더 많은 경향이 있으며, 침실 수보다는 주택의 크기가 주택 가격에 더 큰 영향을 준다. 똑같은 크기의 두 집이 있다고 하면, 작은 크기의 침실이 여러 개 있는 것을 선호하지 않는 것이 합리적이다.
이렇게 상호 연관된 예측변수들을 사용하면 회귀계수의 부호와 값의 의미를 해석하기가 어려울 수 있다. 침실 수, 평수, 욕실 수에 대한 변수들은 모두 상관관계가 있다. 이 세 변수를 제거한 후 모형은
update(step_lm, . ~ . -SqFtTotLiving - SqFtFinBasement - Bathrooms) 침실 수에 대한 계수가 양수로 바뀐다.
상호 연관된 변수들은 회귀 계수를 해석할 때 한가지 문제일 뿐이다. house_lm에는 주택의 위치를 설명하는 변수가 없어서 모델은 서로 다른 유형의 지역이 섞여 있다. Location은 confounding variable일 지도 모른다.
회귀분석에서는 다중공선성 문제를 반드시 해결해야 한다. 다중공선성이 사라질 때까지 변수를 제거해야 한다. perfect multicollinearity는 패키지들이 알아서 해결을 해 주지만, nonperfect multicollinearity인 경우는 결과가 불안정할 수 있다.
다중공선성은 tree, clustering, and nearest-neighbors같은 방법에서는 그다지 문제가 되지 않으며, 이들 방법에서는 P-1개 대신에 P개의 dummy variable을 유지하는 것이 좋다. 물론 이러한 방법에서도 예측변수의 비중복성은 유지하자.
house_lm을 다시 생각해보자. SqFtLot, Bathrooms, and Bedrooms는 모두 음수였다. 이전에 만든 ZipGroup을 포함해보자. 회귀 계수를 보면 ZipGroup이 아주 중요한 변수임을 알 수 있고 SqFtLot and Bathrooms은 이제 양수이며 Bedroom은 여전히 음수이다. 그렇게 직관적이지는 않지만, 이것은 부동산 업계에서 잘 알려진 현상이다. 살기 좋은 지역에서는 욕실 수가 같은 주택의 경우 작은 침실이 여러 개 있으면 오히려 값어치가 떨어진다.
모델에서 주효과만 사용한다면, 여기에는 예측변수와 응답변수 간의 관계가 다른 예측변수들에 대해 독립적이라는 암묵적인 가정이 있다. 하지만 이것은 종종 사실이 아니다.
예를 들어, King County 주택 데이터에서, 부동산에서 위치가 가장 중요하다는 사실은 모두 알고 있다. 주택 크기와 매매 가격 간의 관계가 위치에 달려 있다고 가정하는 것은 자연스러운 일이다. 임대료가 싼 지역에 지어진 큰 집은 비싼 지역에 지어진 큰 집과 같은 가치를 유지하기가 어려울 것이다.
Model Selection with Interactions Terms
주요 개념
이상치는 사기 사건이나 갑작스러운 사건 발생과도 관련이 있다. 따라서 이상치를 발견하는 것은 아주 중요한 비즈니스 니즈가 될 수 있다.
회귀모형에서 제외됐을 때 중요한 변화를 가져오는 값을 영향치라고 한다. 회귀분석에서, 잔차가 크다고 해서 모두 이런 값이 되는 것은 아니다. hat-value가 2(p + 1)/n 이상인 값들은 레버리지가 높은 데이터 값을 나타낸다. Cook’s distance는 레버리지와 잔차의 크기를 합쳐서 영향력을 판단한다. 쿡의 거리가 4 / (n - p - 1)보다 크면 영향을 크게 미친다.
회귀모형을 적합하는 목적이 새로 들어오는 값에 대해 믿을 만한 예측값을 얻기 위함이라면, 데이터의 크기가 작을 경우에만 영향력이 큰 관측 데이터를 확인하는 작업이 유용하다. 데이터수가 클 경우, 어떤 한 값이 회귀방정식에 엄청난 변화를 가져오기란 쉽지 않다. 물론 이상 검출이 목적이라면 높은 영향력의 값들을 찾는 것이 도움이 된다.
통계학자들이 비선형 회귀에대해 이야기할 때 최소제곱법으로 적합할 수 없는 모델을 나타낸다. 설명변수들 또는 변환한 설명변수들의 선형 결합으로 표현될 수 없는 반응변수를 가지는 모델을 말한다. 비선형모델은 수치적 최적화를 필요로 하기 때문에 적합이 어렵고 더 많은 계산을 필요로 한다. 이러한 이유로 일반적으로 가능하면 선형 모형을 선호한다.
A polynomial regression fit for the variable SqFtTotLiving(solid line) versus a smooth (dashed line; see the following section about splines)
spline에 대한 기술적 정의는 일련의 조각별 연속 다항식을 뜻한다. 구간별 다항식은 예측변수를 위한 일련의 고정된 점(knot)사이를 부드럽게 연결한다. spline을 구하는 것은 polynomail regression보다 훨씬 복잡하다. cf. Introduction to spline2
library(splines)
knots <- quantile(house_98105$SqFtTotLiving, p=c(.25, .5, .75))
lm_spline <- lm(AdjSalePrice ~ bs(SqFtTotLiving, knots=knots, degree=3) +
SqFtLot + Bathrooms + Bedrooms + BldgGrade, data=house_98105)다항식의 차수와 knot를 설정해야 한다.
선형 항에서는 계수가 변수에 대한 직접적인 의미를 갖는 반면, spline 항에서는 해석이 힘들다. 대신, spline의 적합도를 확인하기 위해 시각화를 하자. polynomial에 비해 spline이 좀 더 매끄럽게 매칭되며 유연성이 뛰어나다. 하지만 이것이 spline regression 이 더 좋다는 것을 의미하지는 않는다. 크기가 아주 작은 주택이 약간 큰 주택보다 더 높은 가치를 가질 것이라는 예측 결과는 경제적으로 맞지 않다. 이것은 confounding variable 때문일 수 있다.
A spline regression fit for the variable SqFtTotLiving (solid line) compared to a smooth (dashed line)
GAM은 spline회귀를 자동으로 적합하는 기술이다.
library(mgcv)
lm_gam <- gam(AdjSalePrice ~ s(SqFtTotLiving) + SqFtLot + Bathrooms + Bedrooms + BldgGrade,
data=house_98105)A GAM regression fit for the variable SqFtTotLiving (solid line) compared to a smooth (dashed line)
주요 개념
More than two categories?
결과가 세 개 이상인 경우도, 조건부확률을 사요하여 여러 개의 이진 문제로 돌려서 생각해볼 수 있다. 예를 들면 앞의 재계약의 문제의 경우도, 두 가지 이진 예측 문제로 볼 수 있다.
이 경우 문제를 두 가지 경우, 고객이 계약을 해지하는 경우와 재계약에 동의하는 경우(동의할 경우 어떤 재계약을 원하는지)로 나누는 것이 좋다. 모델 피팅 관점에서도, 다중 분류 문제를 일련의 이진 문제로 변환하는 것이 종종 유리하다. 하나의 범주가 다른 범주보다 훨씬 더 일반적인 경우 특히 그렇다.
is_loan_char <- vapply(loan_data, is.character, logical(1))
loan_data[, is_loan_char] <- lapply(loan_data[, is_loan_char], as.factor)
library(klaR)## Warning: package 'klaR' was built under R version 3.5.2
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
naive_model <- NaiveBayes(outcome ~ purpose_ + home_ + emp_len_, data = na.omit(loan_data))
naive_model$tables # P(X_j | Y = i)## $purpose_
## var
## grouping credit_card debt_consolidation home_improvement major_purchase
## default 0.15151515 0.57571347 0.05981209 0.03727229
## paid off 0.18759649 0.55215915 0.07150104 0.05359270
## var
## grouping medical other small_business
## default 0.01433549 0.11561025 0.04574126
## paid off 0.01424728 0.09990737 0.02099599
##
## $home_
## var
## grouping MORTGAGE OWN RENT
## default 0.4313440 0.0832782 0.4853778
## paid off 0.4894800 0.0808963 0.4296237
##
## $emp_len_
## var
## grouping < 1 Year > 1 Year
## default 0.04728508 0.95271492
## paid off 0.03105289 0.96894711
naive Bayesian classifier은 biased 추정치를 구한다. 하지만, 목적이 Y = 1일 확률에 따라 레코드들에 순위를 메기는 것이 목적이므로 확률의 비편향된 추정치가 필요로하지 않아서 나이브베이즈는 나름 우수한 결과를 보인다.
# MASS는 dplyr의 select와 충돌
loan_lda <- MASS::lda(outcome ~ borrower_score + payment_inc_ratio, data=loan3000)
loan_lda$scaling## LD1
## borrower_score 7.17583880
## payment_inc_ratio -0.09967559
pred <- predict(loan_lda)
head(pred$posterior)## default paid off
## 1 0.5535437 0.4464563
## 2 0.5589534 0.4410466
## 3 0.2726962 0.7273038
## 4 0.5062538 0.4937462
## 5 0.6099525 0.3900475
## 6 0.4107406 0.5892594
lda_df <- cbind(loan3000, prob_default=pred$posterior[,'default'])
ggplot(data=lda_df, aes(x=borrower_score, y=payment_inc_ratio, color=prob_default)) +
geom_point(alpha=.6) +
scale_color_gradient2(low='white', high='blue') +
geom_line(data=lda_df0, col='green', size=2, alpha=.8) 일부 응용 분야에서 LAD를 사용하고 있으며, PCA와 같이 아직도 많이 사용되는 다른 방법들과도 연결된다. 설명변수의 중요도나 feature selectiion에도 사용된다.
다변량 분석에서 변수들을 표준화하기 위해 공분산 행렬 사용. 마할라노비스 거리(Mahalanobis distance)라고 부른다.
피셔의 선형판별: \(SS_{between} / SS_{within}\)을 최대화하는 설명변수들의 선형결합을 찾는다.
희귀한 클래스에 속한 구성원을 확인하는 것이 목표라면 이 기준값을 더 낮게 할수록 좋은 경우가 종종 있다.
선형회귀에 적용되었던 많은 개념이 로지스틱 회귀(그리고 다른 GLM들)에도 똑같이 이어진다. 예를 들면 여기에도 stepwise, interaction term, spline term등을 모두 사용할 수 있다. confounding variable이나 correlated variable과 관련한 문제들도 동일하게 고려해야 한다.
# GAM도 사용 가능
logistic_gam <- gam(outcome ~ s(payment_inc_ratio) + purpose_ + home_ + emp_len_ + s(borrower_score),
data=loan_data, family='binomial')
# partial residual plot
terms <- predict(logistic_gam, type='terms')
partial_resid <- resid(logistic_model) + terms
df <- data.frame(payment_inc_ratio = loan_data[, 'payment_inc_ratio'],
terms = terms[, 's(payment_inc_ratio)'],
partial_resid = partial_resid[, 's(payment_inc_ratio)'])
ggplot(df, aes(x = payment_inc_ratio, y = partial_resid, solid = FALSE)) +
geom_point(shape = 46, alpha=.4) +
geom_line(aes(x = payment_inc_ratio, y = terms),
color = 'red', alpha = .5, size = 1.5) +
labs(y = 'Partial Residual')Partial residuals from logistic regression
위쪽 구름은 1의 응답(연체)을 의미하고 아래쪽 구름은 0의 응답(대출 상환)을 의미한다. 결과변수가 이진형이기 때문에, 로지스틱 회귀에서 얻은 잔차는 보통 이러한 형태를 띄게 된다. 로지스틱 회귀에서 편잔차는 회귀에서보다 덜 중요하긴 하지만, 비선형성을 검증하고 영향력이 큰 레코드들을 확인하는데 여전히 유용하다.
pred <- predict(logistic_gam, newdata=train_set)
pred_y <- as.numeric(pred > 0)
true_y <- as.numeric(train_set$outcome=='default')
true_pos <- (true_y==1) & (pred_y==1)
true_neg <- (true_y==0) & (pred_y==0)
false_pos <- (true_y==0) & (pred_y==1)
false_neg <- (true_y==1) & (pred_y==0)
conf_mat <- matrix(c(sum(true_pos), sum(false_pos),
sum(false_neg), sum(true_neg)), 2, 2)
colnames(conf_mat) <- c('Yhat = 1', 'Yhat = 0')
rownames(conf_mat) <- c('Y = 1', 'Y = 0')
conf_matConfusion matrix for a binary response and various metrics
cf. F1-Score
cf. 분류모델의 성능 평가
일반적으로 1이 더 중요한 사건이다. 예를 들면 사기성 보험 청구를 정확히 잡아내는 것은 몇 천 달러의 돈을 아끼는 결과를 가져다준다. 반대로 사기성이 아닌 보험 청구를 정확하게 파악하는 것은, 사기성으로 의심되는 보험 청구를 일일이 손으로 확인하는 데 드는 비용과 노력만을 절약해줄 뿐이다.
클래스를 쉽게 분리하기 어려운 경우에는, 가장 정확도가 높은 모델은 모든 것을 무조건 0으로 분류하는 모델일 수도 있다. 그러나 이 모델은 결국 있으나 마나다. 대신, 전반적인 정확도가 떨어지더라도 1을 잘 분류해 내는 모델이 있다면 그 모델을 선호할 것이다.
단순히 cutoff를 희귀 클래스의 비율로 정할 수도 있다.
Receiver Operation Characteristics
idx <- order(-pred)
recall <- cumsum(true_y[idx] == 1) / sum(true_y == 1)
specificity <- (sum(true_y == 0) - cumsum(true_y[idx] == 0)) / sum(true_y == 0)
roc_df <- data.frame(recall = recall, specificity = sepcificity)
ggplot(roc_df, aes(specificity, recall)) +
geom_line(color = "blue") +
scale_x_reverse(expaned = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
geom_line(data = data.frame(x = (0:100) / 100), aes(x, 1 - x), linetype = "dotted", color = "red")클래스 불균형이 심하다면 ROC curve대신에 PR(precision-recall) curve
cf. Business Value
cf. Lift Analysis
모델의 전반적인 성능을 알 수 있다. slope of lift chart가 not monotonic이면 결함을 빠르게 찾을 수 있다.
분류기 성능을 평가하는 지표로 AUC를 사용하면 단순히 정확도만을 사용하는 것보다는 나은 결과를 얻을 수 있따. 전체적인 정확도도 높이면서 실무에서 중요한 1을 더 정확히 분류해야 하는 trade-off를 얼마나 잘 처리하는지 평가할 수 있기 때문이다. 하지만 희귀 케이스 문제에서는 모든 레코드를 0으로 분류하지 않도록 하려면 모델의 확률 cut-off를 0.5미만으로 낮춰야 하는 문제가 있따. 즉, 1의 중요성을 너무 크게 반영하여 1을 과대평가하는 결과를 낳을 수 있다.
리프트 차트를 계산하려면 y축에 recall, x축에 총 레코드 수를 나타내는 cumulative gains chart를 작성해야 한다.
리프트 곡선은 레코드를 1로 분류하기 위한 확률 컷오프 값에 따른 결과의 변화를 한눈에 볼 수 있게 해준다.
uplift
가끔 리프트와 동일한 의미로 uplift라는 표현을 사용한다. 하지만 이 용어는 좀 더 제한적인 상황에서 사용된다. A/B검정을 수행하고, 처리 A나 B 가운데 하나를 예측변수로 사용하는 예측 모델에서 사용된다. 업리프트는 처리 A와 처리 B 사이의 개별적인 한 케이스에 대해 예측된 결과의 향상을 의미한다. 처음에는 예측변수를 A로 놓고 개별 데이터를 점수화하고 그리고 다시 예측 변수를 B로 바꿔서 한 다음 점수를 보고 결정한다. 영업 담당자나 정치 선거 컨설턴트가 고객이나 유권자를 대상으로 한 두 메시지 중에 어느 것이 더 효과적인지를 결정하는 데 등에 이 방법을 사용한다.
다수의 클래스에 해당하는 데이터에 대한 정보를 더 많이 사용하여 소수 클래스에 대한 예측 값이 과소추정된다. 일반적으로 소수 클래스의 데이터가 수만 개 정도 있따면 충분하다. 물론 1과 0을 분리하기가 쉽다면, 더 적은 데이터로 충분할 수도 있다.
undersampling의 단점으로 일부 데이터가 버려지기 떄문에 모든 정보를 활용하지 못한다는 점이다. 상대적으로 작은 데이터 집합에서, 희귀 클래스 경우의 레코드가 몇백 혹은 몇천 개라면, 다수 클래스에 대한 과소표본추출은 정말 유용한 정보까지 버리게 되는 결과를 초래할 수 있다. 이럴 경우, 복원추출 방식(bootstrap)으로 희귀 클래스의 데이터를 Oversampling해야 한다.
데이터에 가중치를 적용하는 방식으로 이와 비슷한 효과를 얻을 수 있다. 많은 분류 알고리즘에서 상향/하향 가중치를 데이터에 적용하기 위해 weight라는 인수를 지원한다. 예를 들면 glm함수에서 weight라는 인수를 사용해서 대출 데이터에 가중치 벡터를 적용해보자.
wt <- ifelse(loan_all_data$outcome =='default',
1/mean(loan_all_data$outcome == 'default'), 1)
full_model <- glm(outcome ~ payment_inc_ratio + purpose_ + home_ + emp_len_+ dti + revol_bal + revol_util, data=loan_all_data, weight=wt, family='binomial')
pred <- predict(full_model)
mean(pred > 0)연체에 대한 가중치를 1/p로 두었다. 여기서 p는 연체의 확률값이다. 그리고 대출 상환에 대한 가중치는 1로 두었다. 연체와 상환의 가중치 합은 거의 동일하다.
가중치를 적용하는 방식은 up/under sampling 모두 적용할 수 있다.
Adapting the Loss Function
많은 분류/회귀 알고리즘은 어떤 기준 혹은 loss function을 최적화한다. 예를 들어, 로지스틱 회귀분석은 deviance를 최소화한다. 어떤 자료에서는 희귀 클래스에 의해 발생하는 문제를 피하기 위한 loss function의 수정을 제안한다. 실무에서는 이것이 어렵다. 분류 알고리즘은 복잡하고 어려워서 수정하기 어렵다. 가중치를 두는 것은 가중치가 높은 데이터를 선호하고 가중치가 낮은 데이터의 오류를 줄여주는 식으로 loss function를 변경하는 쉬운 방법이다.
부트스트랩을 통한 업샘플링 방식의 변형으로 기존에 존재하는 데이터를 살짝 바꿔 새로운 데이터를 만드는 data generation 방법이 있다. 이 방법에는 데이ㅌ의 개수가 제한적이기 떄문에 알고리즘을 통해 분류 “규칙”을 세우기에는 정보가 충분하지 않다는 직관이 바탕에 깔려 있다. 비슷하지만 기존의 데이터와 다른 데이터를 생성해서 좀 더 로버스트한 분류 규칙을 배울 수 있는 기회를 주고자 하는 것이다. 이는 앙상블 모델에 담겨 있는 개념과 매우 비슷하다.
SMOTE(Synthetic Minority Oversampling Technique)알고리즘은 업샘플링된 레코드와 비슷한 레코드를 찾고, 원래 레코드와 이웃 레코드의 랜덤 가중평균으로 새로운 합성 레코드를 만든다. 여기에 대해 각각의 예측변수에 대해 개별적으로 가중치를 생성한다. 새로 합성된 업샘플 레코드의 개수는 데이터의 균형을 맞추기 위해 필요한 업샘플링 비율에 따라 달라진다.
신규 대출에서 연체로 인해 예상되는 비용이 C라고 하고 대출 상환을 통해 얻을 수 있는 수익을 R라고 하자. 이때 신규 대출의 기대 수익은
\(기대 수익 = P(Y = 0)*R + P(Y = 1)*C\)
여기서 대출 결과를 단순히 연체나 상환, 둘 중 하나로 결정하는 대신에, 대출을 통해 얻을 수 있는 기대 수익이 있는지 없는지를 결정하는 것이 더 말이 된다. 대출을 갚지 않을 확률을 예측하는 것은 중간 단계이고, 결국은 사업의 목적인 기대 수익을 결정하기 위해 대출의 전체적인 가치를 얻어내야 한다. 예를 들면, 가치가 적은 대출보다는 연체 확률이 더 높더라도 가치가 더 큰 대출을 선호하는 편이 나을 수도 있는 것이다.
기억: 예측값에 대한 탐색적 분석은 필요하다.
Distance Metrics
One Hot Encoder
이진 가변수의 집합으로 변환한다. 원핫인코딩이라는 용어는 디지털 회로 분야에서 유래한 것으로, 하나의 비트만 양수(hot)가 허용되는 회로 설정을 말한다.
loan_std <- scale(loan_df)
knn_pred <- knn(train=loan_std, test=newloan_std, cl=outcome, k=5)
loan_df[attr(knn_pred,"nn.index"),]각 변수가 예측력 측면에서 갖는 중요성이 모두 같다는 것을 의미한다. 일부 변수가 다른 변수보다 중요하다는 주관지 지식이 있는 경우, 조정할 수 있따. 예를 들면 대출 데이터의 경우, 소득 대비 상환 비율이 매우 중요하다고 판단하는 것이 합리적이다.
일반적으로 K가 너무 작으면 데이터의 노이즈 성분까지 고려하는 오버피팅 문제가 발생한다. K값이 클수록 훈련 데이터에서 오버피팅의 리스크를 감소시키는 smoothing 제공한다. 반면에, K가 너무 크면 데이터를 oversmooth해서 주요 이점중 하나인 데이터의 지역 정보를 포착하는 능력을 놓친다.
CV를 사용하여 K를 결정하자. 데이터에 노이즈가 거의 없고 아주 잘 구조화된 데이터의 경우 K값이 작을수록 잘 동작한다. 신호처리 분야에서 사용하는 전문 용어를 잠시 빌려오자면, signal-to-noise ratio(SNR)가 높은 데이터라는 의미이다. 일반적으로 손글씨 데이터 또는 음성 데이터가 SNR이 높은 데이터라고 할 수 있다. 반면에 대출 데이터와 같이 노이즈가 많아 SNR이 낮은 데이터의 경우, K가 클수록 좋다. 보통 K를 1에서 20사이에 놓는다. 동률이 나오는 경우를 막기 위해 보통은 홀수를 사용한다.
Bias-Variance Tradeoff
분산은 학습 데이터 선정에 따라 발생하는 모델링 오차를 말한다. 즉 다른 학습 데이터를 사용할 경우 모델이 달라지는 정도를 의미한다. 편향은 실제 세계를 잘 반영하지 못하기 떄문에 일어나는 모델링 오차를 의미한다. 이러한 오차는 추가적인 훈련 데이터셋을 더해도 사라지지 않는다. 유연한 모델이 오버피팅 될 때 분산은 증가한다. CV를 사용하자.
다른 분류 방법들의 특정 단계에서 사용할 수 있게 local knowledge를 추가하기 위해 KNN을 사용할 수 있다.
여기서 다중공선성은 문제가 되지 않는다. 위의 2단계에서 얻은 정보는 소수의 근접한 레코드들로부터 얻은 매우 지협적인 정보이기 때문이다. 따라서 새로 얻은 정보는 불필요하거나, 중복석이 있지 않다.
이러한 단계의 KNN사용은 다수의 예측 모델 방법은 서로서로 함께 사용되는 앙상블 학습의 형태로 볼 수 있다.
borrow_df <- model.matrix(~ -1 + dti + revol_bal + revol_util + open_acc + delinq_2yrs_zero + pub_rec_zero, data=loan_data)
borrow_knn <- knn(borrow_df, test=borrow_df, cl=loan_data[, 'outcome'], prob=TRUE, k=10)
prob <- attr(borrow_knn, "prob")
borrow_feature <- ifelse(borrow_knn=='default', prob, 1-prob)
summary(borrow_feature)용어 정리
트리 모델이란 if-then-else 규칙의 집합체라 할 수 있다. 단순 회귀나 로지스틱 회귀와는 반대로 트리는 데이터에 복잡한 상호작용에 대응하는 숨겨진 패턴을 발견하는 능력이 있다. KNN이나 naive Bayes와는 달리 단순한 트레 모델은 쉽게 해석 가능한 예측변수의 관계로 표현될 수 있다.
cf. https://cran.r-project.org/web/packages/rpart/vignettes/longintro.pdf
library(rpart)
loan_tree <- rpart(outcome ~ borrower_score + payment_inc_ratio, data = loan_data, control = rpart.control(cp = .005)) # cp의 default는 0.01
plot(loan_tree, uniform = TRUE, .05)
text(loan_tree)
loan_treeThe rules for a simple tree model fit to the loan data
트리의 깊이를 들여쓰기의 정도로 쉽게 파악할 수 있다. 각 노드는 해당 분할 규칙에 대해 우세한 쪽으로 결정하는 임시 분류를 의미한다. ’손실’은 이러한 임시 분할에서 발생하는 오분류의 개수를 의미한다. 노드, (실제)분류, (실제)오분류, (예측한)분류, 해당 데이터에서 클래스별 비율.
의사 결정 트리를 만들 때는 재귀 분할이라고 하는 알고리즘을 사용한다. 예측변수 값을 기준으로 데이터를 반복적으로 분할해나간다.
The rules for a simple tree model fit to the loan data
알고리즘의 재귀 부분은 1. 전체 데이터를 가지고 A를 초기화한다. 2. A를 두 부분 \(A_1\)과 \(A_2\)로 나누기 위해 분할 알고리즘을 적용한다. 3. \(A_1\)과 \(A_2\) 각각에서 2번 과정을 반복한다. 4. 분할을 해도 더 이상 하위 분할 영역의 동질성이 개선되지 않을 정도로 충분히 분할을 진행했을 때, 알고리즘을 종료한다.
\(P(Y = 1) = \frac{Number of 1s in the partition}{Size of the partition}\)
accuracy는 impurity를 측정하는 데에 썩 좋지 않으므로 Gini impurity.와 entropy가 대표적인 impurity 측정 지표이다. 클래스가 2개 이상인 분류 문제에도 적용가능하다.
주의: gini impurity와 gini coefficient를 혼동하지 말자. 지니 계수는 이진 분류 문제로 한정되며, AUC지표와 관련이 있는 용어다.
앞에서 설명한 분할 알고리즘에서 이 불순도 측정 지표를 사요했다. 분할로 만들어지는 각 영역에 대해 불순도를 측정한 다음, 가중평균을 계산하고 (단계마다) 가장 낮은 가중평균을 보이는 분할 영역을 선택했다.
Pruning
holdout data에서의 에러가 최소가 되는 지점까지 가지치기를 진행한다. 여러 트리에서 얻은 예측 결과를 하나로 합칠 때에도 마찬가지로 성장을 멈추는 방법이 필요하다. 가지치기는 앙상블 기법에서 사용할 트리들을 얼마나 크게 만들지 결정하기 위한 교차타당성검사에서 중요한 역할을 담당한다.
가지 분할을 멈추는 대표적인 두 가지 방법
다소 임의적이라고 할 수 있는 첫 번째 방법은 탐색 작업에 유용할 수 있지만 최적값을 결정하기가 매우 어렵다. cp를 이용하면 어떤 크기의 트리가 새로운 데이터에 대해 가장 좋은 성능을 보일지 추정할 수 있다.
cp가 매우 작다면 트리는 실제 의미 있는 신호뿐 아니라 노이즈까지 학습하여 오버피팅되는 문제가 발생하게 될 것이다. 반면에 cp가 너무 크다면 트리가 너무 작아 예측 능력을 거의 갖지 못할 것이다. rpart함수의 기본 설정은 0.01로 다소 큰 값이다. 이전 예제에서 기본값을 그대로 사용하면 트리 분할이 단 한 번만 일어나므로, 예제에서는 cp를 0.005로 지정했다. 보통은 사전 탐색 분석에서 몇 가지 값을 테스트해보는 것으로도 충분하다.
최적의 cp를 결정하는 것은 bias-variance tradeoff를 보여주는 하나의 대표적인 예라고 할 수 있다. cp를 추정하는 가장 일반적인 방법은 교차타당성 검정을 이용하는 것이다.
차이점은 각 하위 분할 영역에서 squared deviations from the mean(squared errors) 값을 이용해 불순도를 측정한다는 점과 RMSE를 이용해 예측 성능을 평가한다는 점이다.
트리 모델은 아래 두 가지 장점을 갖는다.
앙상블 방법의 간단한 버전 1. 주어진 데이터에 대해 예측 모델을 만들고 예측 결과를 기록한다. 2. 같은 데이터에 대해 여러 모델을 만들고 결과를 기록한다. 3. 각 레코드에 대해 예측된 결과들의 평균(또는 가중평균, 다수결 투표)을 구한다.
Bootstrap aggregating
랜덤 포레스트는 의사 결정 나무 모형에서 한 가지 중요한 요소가 추가된 배깅 방법을 적용한 모델이다. 레코드를 표본추출할 때, 변수 역시 샘플링하는 것이다. 일반적인 의사 결정 트리에서는 하위 분할 영역 A를 만들 때, 지니 불순도와 같은 기준값이 최소화되도록 변수와 분할지점을 결정했다. 랜덤 포레스트에서는 알고리즘의 각 단계마다, 고를 수 있는 변수가 랜덤하게 결정된 전체 변수들의 부분집합에 한정된다.
각 단계에서 샘플링하는 변수의 수는 보통 \(\sqrt P\).
library(randomForest)## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
##
## Attaching package: 'randomForest'
## The following object is masked from 'package:dplyr':
##
## combine
## The following object is masked from 'package:ggplot2':
##
## margin
loan3000$outcome <- as.factor(loan3000$outcome)
(rf <- randomForest(outcome ~ borrower_score + payment_inc_ratio, data=loan3000))##
## Call:
## randomForest(formula = outcome ~ borrower_score + payment_inc_ratio, data = loan3000)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 1
##
## OOB estimate of error rate: 39.1%
## Confusion matrix:
## default paid off class.error
## default 869 576 0.3986159
## paid off 597 958 0.3839228
여기서는 예측 변수를 두 개만 사용했기 때문에, 알고리즘은 매 단계마다 둘 중 한 변수를 랜덤하게 선택하게 된다.
The out-of-bag(OOB) estimate of error is the error rate for the trained models, applied to the data left out of the training set for that tree.
error_df <- tibble(error_rate = rf$err.rate[, "OOB"], num_trees = 1:rf$ntree)
ggplot(error_df, aes(num_trees, error_rate)) +
geom_line() + theme_bw()
Figure 2: The improvement in accuracy of the random forest with the addition of more trees
pred <- predict(loan_lda)
pred_df <-
pred$posterior %>%
as.tibble() %>%
mutate(class = factor(pred$class, level = c("paid off", "default")))
rf_df <- cbind(loan3000, pred_df)
ggplot(rf_df, aes(borrower_score, payment_inc_ratio, color = class, shape = class)) +
geom_point(alpha = .6, size = 2) +
scale_shape_manual(values = c(46, 4)) +
theme_bw()
Figure 3: The predicted outcomes from the random forest applied to the loan default data
일반적이지 않은 예외 사항까지 학습해서 생기는 결과로, 랜덤 포레스트에 의한 오버피팅의 위험성을 보여준다.
rf_all <- randomForest(outcome ~ . -status -home_ownership -purpose -emp_length , data = loan_data, importance = TRUE)
rf_all
varImpPlot(rf_all, type=1)
varImpPlot(rf_all, type=2)gini impurity를 굳이 고려할 필요가 있을까? 기본적으로 randomforest함수는 지니 불순도만을 활용한다. 모델의 정확도는 추가적인 계산(랜덤 순열 조합으로 데이터를 뽑고 이 데이터를 이용해 예측 결과를 구하는 일련의 과정)이 필요한 반면, 지니 불순도는 알고리즘상에서 부차적으로 얻어지는 결과물이다. 예를 들어 몇천 개의 모델을 만들어야 해서 계산 복잡도가 중요하다면, 프로덕션 상황에서는 계산을 더 해서 얻는 이득이 거의 없을 수 있다. 또한 지니 불순도의 감소를 통해 분할 규칙을 만드는 데 어떤 변수를 사용하는 것이 좋을지를 더욱 분명히 알 수 있다. 물론 간단한 트리 모델이 갖고 있는 쉽게 해석이 가능하다는 장점을 랜덤 포레스트에서는 찾아볼 수 없다. 지니 계수 감소와 모델 정확도 감소 사이의 차이를 잘 조사해보면, 변수 중요도를 통해 모델을 개선할 수 있는 몇 가지 아이디어들을 얻을 수 있을 것이다.
오버피팅을 피하기 위해서 매우 중요하다
랜덤 포레스트를 노이즈가 많은 데이터에 적용할 때, 기본 설정으로는 오버피팅에 빠질 수 있다. nodesize와 maxnodes를 크게 하면 더 작은 트리를 얻게 되고 거짓 예측 규칙들을 만드는 것을 피할 수 있게 된다. 하이퍼파라미터에 다른 값들을 적용했을 때의 효과를 알아보려면 CV를 이용할 수 있다.
용어 정리
파라미터를 잘 선택하면 랜덤 포레스트를 그대로 모방(emulate) 할 수 있다.
Adaboost의 알고리즘
The boosted estimate is given by: \[\hat F = \sum_{i = 1}^M\alpha_i\hat f_i\]
잘못 분류된 관측 데이터에 가중치를 증가시킴으로써, 현재 성능이 제일 떨어지는 데이터에 대해 더 집중해서 학습을 하도록 하는 효과를 가져온다. \(\alpha_m\)값을 이용해 모델의 오차가 낮을수록 더 큰 가중치를 부여한다.
gradient boosting은 adaboosting과 거의 비슷하지만, 비용함수를 최적화하는 접근법을 사용했다는 점에서 차이가 있다. gradient boosting에서는 가중치를 조정하는 대신에 모델이 pseudo-residual를 학습하도록 한다. 이는 잔차가 큰 데이터를 더 집중적으로 학습하는 효과를 가져온다. stochastic gradient boosting에서는 랜덤 포레스트에서와 유사하게, 매 단계마다 데이터와 예측변수를 샘플링하는 식으로 그레디언트 부스팅에 랜덤한 요소를 추가한다.
library(xgboost)##
## Attaching package: 'xgboost'
## The following object is masked from 'package:dplyr':
##
## slice
loan_data$X1 <- NULL
predictors <- data.matrix(loan3000[, c('borrower_score', 'payment_inc_ratio')])
label <- as.numeric(factor(loan3000$outcome, level = c("paid off", "default"))) - 1
xgb <- xgboost(data=predictors, label=label,
objective = "binary:logistic",
params=list(subsample=.63, eta=0.1), nrounds=100)## [1] train-error:0.366333
## [2] train-error:0.352667
## [3] train-error:0.342667
## [4] train-error:0.338000
## [5] train-error:0.336000
## [6] train-error:0.333000
## [7] train-error:0.329667
## [8] train-error:0.324667
## [9] train-error:0.323667
## [10] train-error:0.322667
## [11] train-error:0.318000
## [12] train-error:0.315000
## [13] train-error:0.315667
## [14] train-error:0.312667
## [15] train-error:0.312333
## [16] train-error:0.313667
## [17] train-error:0.312000
## [18] train-error:0.314000
## [19] train-error:0.315000
## [20] train-error:0.315000
## [21] train-error:0.311000
## [22] train-error:0.309333
## [23] train-error:0.304000
## [24] train-error:0.303667
## [25] train-error:0.301667
## [26] train-error:0.303333
## [27] train-error:0.303000
## [28] train-error:0.304000
## [29] train-error:0.299333
## [30] train-error:0.301667
## [31] train-error:0.300667
## [32] train-error:0.299000
## [33] train-error:0.300000
## [34] train-error:0.296333
## [35] train-error:0.298333
## [36] train-error:0.297333
## [37] train-error:0.297000
## [38] train-error:0.299000
## [39] train-error:0.294000
## [40] train-error:0.292000
## [41] train-error:0.293667
## [42] train-error:0.292333
## [43] train-error:0.286667
## [44] train-error:0.286000
## [45] train-error:0.282667
## [46] train-error:0.280667
## [47] train-error:0.281000
## [48] train-error:0.278000
## [49] train-error:0.276667
## [50] train-error:0.274333
## [51] train-error:0.274667
## [52] train-error:0.274000
## [53] train-error:0.274000
## [54] train-error:0.276000
## [55] train-error:0.275000
## [56] train-error:0.273667
## [57] train-error:0.272667
## [58] train-error:0.272000
## [59] train-error:0.272000
## [60] train-error:0.272000
## [61] train-error:0.273000
## [62] train-error:0.270667
## [63] train-error:0.271333
## [64] train-error:0.268333
## [65] train-error:0.265667
## [66] train-error:0.265667
## [67] train-error:0.265000
## [68] train-error:0.264667
## [69] train-error:0.266000
## [70] train-error:0.264667
## [71] train-error:0.260333
## [72] train-error:0.261667
## [73] train-error:0.259000
## [74] train-error:0.256000
## [75] train-error:0.253667
## [76] train-error:0.252000
## [77] train-error:0.250667
## [78] train-error:0.249000
## [79] train-error:0.248333
## [80] train-error:0.248333
## [81] train-error:0.246667
## [82] train-error:0.247667
## [83] train-error:0.246000
## [84] train-error:0.245667
## [85] train-error:0.245667
## [86] train-error:0.246333
## [87] train-error:0.245333
## [88] train-error:0.247000
## [89] train-error:0.247667
## [90] train-error:0.245667
## [91] train-error:0.245000
## [92] train-error:0.238000
## [93] train-error:0.237667
## [94] train-error:0.237667
## [95] train-error:0.237000
## [96] train-error:0.239000
## [97] train-error:0.237667
## [98] train-error:0.238000
## [99] train-error:0.237000
## [100] train-error:0.237333
pred <- predict(xgb, newdata = predictors)
xgb_df <- cbind(loan3000, pred_default = pred >.5, prob_default = pred)
ggplot(xgb_df, aes(borrower_score, payment_inc_ratio, color = pred_default, shape = pred_default)) +
geom_point(alpha = .6, size = 2)
Figure 4: The predicted outcomes from XGBOost applied to the loan default data
전체적으로 랜덤 포레스트를 이용했을 때와 예측 결과가 비슷한 것을 볼 수 있다. 예측 결과가 다소 복잡하고, 아직도 신용점수가 높은 경우인데도 가끔 예측 결과가 연체로 나오기도 한다.
Blind application of xgboost can lead to unstable models as a result of overfitting to the training data. The problem with overfitting is twofold:
어떤 모델링 기법이든 오버피팅에 바지기 쉽다. 예를 들어 회귀에서 너무 많은 변수를 고려하다 보면 결국 모델은 잘못된 예측 결과를 만들어 낼 수 있다. 하지만 대부분의 통계적인 기법들은 예측변수들을 까다롭게 걸러내는 과정을 통해 오버피팅을 피할 수 있다. 심지어 랜덤포레스트는 별다른 파라미터 튜닝 없이도 상당히 괜찮은 모델을 만들곤 한다. 하지만 xgboost의 경우는 그렇지 않다. 대출 데이터에서 모든 변수를 고려해 xgboost를 학습해 보자.
library(Matrix)##
## Attaching package: 'Matrix'
## The following object is masked from 'package:tidyr':
##
## expand
library(caret)## Loading required package: lattice
##
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
##
## lift
#status는 빼자
ohe_feats <- colnames(loan_data[, vapply(loan_data, is.factor, logical(1))])[-c(1, 5)]
dummies <- dummyVars(~., data = loan_data[, ohe_feats])
df_all_ohe <- as.data.frame(predict(dummies, newdata = loan_data))
df_all_combined <-
loan_data %>%
dplyr::select(-which(colnames(loan_data) %in% ohe_feats), -status, -outcome) %>%
bind_cols(df_all_ohe)
predictors <- df_all_combined %>% as.matrix()
label <- as.numeric(factor(loan_data$outcome, level = c("paid off", "default"))) - 1
test_idx <- sample(nrow(predictors), 10000)
xgb_default <- xgboost(data = predictors[-test_idx, ], label=label[-test_idx],
objective = "binary:logistic", nrounds=250, verbose = 0)
pred_default <- predict(xgb_default, predictors[test_idx,])
error_default <- abs(label[test_idx] - pred_default) > 0.5
xgb_default$evaluation_log[250,]## iter train_error
## 1: 250 0.136778
mean(error_default)## [1] 0.3513
전체 데이터 가운데 랜덤하게 테스트를 위한 10000개를 데이터를 선별하고, 나머지 데이터를 학습데이터로 사용했다. 부스팅 모델을 학습한 결과, 학습 데이터에 대한 오차율은 250, 0.136778 였지만, 테스트 데이터에 대한 오차율은 그것보다 훨씬 높은 0.3513 였다. 이는 오버피팅의 결과이다.
앞서 살펴본 것처럼 subsample이나 eta같은 파라미터를 통해 부스팅 시 오버피팅을 방지할 수 있다. 여기에서 소개할 또 다른 방법은 refularization방법이다. 이는 모델의 복잡도에 따라 벌점을 추가하는 형태로 비용함수를 변경하는 방법이다. 의사 결정 트리에서는 지니 불순도와 같은 비용 기준값을 최소화하는 쪽으로 모델을 피팅했다. xgboost에서는 모델의 복잡도를 의미하는 항을 추가하는 형태로 비용함수를 변경할 수 있다.
xgboost에는 모델을 정규화하기 위한 두 파라미터 alpha와 lambda가 존재한다. 이는 각각 맨하탄 거리와 유클리드 거리를 의미한다. 이 파라미터들을 크게 하면, 모델이 복잡해질수록 더 많은 벌점을 부여하게 되고 결과적으로 얻어지는 트리의 크기가 작아지게 된다. 예를 들어 lambda의 값을 1000으로 하면 무슨 일이 벌어지는지 알아보자.
#verbose = 0 이면 message출력 x
xgb_penalty <- xgboost(data=predictors[-test_idx,],
label=label[-test_idx],
params=list(eta=.1, subsample=.63, lambda=1000),
objective = "binary:logistic", nrounds=250, verbose=0)
pred_penalty <- predict(xgb_penalty, predictors[test_idx,])
error_penalty <- abs(label[test_idx] - pred_penalty) > 0.5
xgb_penalty$evaluation_log[250,]## iter train_error
## 1: 250 0.313253
predict함수는 좀 더 편리한 파라미터 ntreelimit를 제공한다. 이는 예측을 위해 첫 i개의 트리 모델만을 사용하는 것을 가능하게 한다. 이를 통해 예측을 위해 사용하는 모델의 개수에 따른 표본 내 오차율과 표본 밖 오차율을 더 쉽게 비교할 수 있다.
error_default <- rep(0, 250)
error_penalty <- rep(0, 250)
for(i in 1:250){
pred_def <- predict(xgb_default, predictors[test_idx,], ntreelimit=i)
error_default[i] <- mean(abs(label[test_idx] - pred_def) >= 0.5)
pred_pen <- predict(xgb_penalty, predictors[test_idx,], ntreelimit = i)
error_penalty[i] <- mean(abs(label[test_idx] - pred_pen) >= 0.5)
}모델 출력에서 xgb_default$evaluation_log은 학습 데이터에 대한 오차를 알려준다. 이를 표본 밖 오차와 결합하여 반복 횟수에 따른 오차율에 대한 그래프를 구할 수 있다.
errors <- bind_rows(xgb_default$evaluation_log,
xgb_penalty$evaluation_log,
data.frame(iter=1:250, train_error=error_default),
data.frame(iter=1:250, train_error=error_penalty))
errors$type <- rep(c('default train', 'penalty train',
'default test', 'penalty test'), rep(250, 4))
ggplot(errors, aes(x=iter, y=train_error, group=type)) +
geom_line(aes(linetype=type, color=type)) +
theme_bw() +
labs(x = "Iterations", y = "Error")
Figure 5: The error rate of the default XGBoost versus a penalized version of XGBoost
RIDGE REGRESSION AND THE LASSO
ridge regression에서는 잔차제곱합에 회귀계수의 개수와 크기에 따라 벌점을 추가한 값을 최소화한다.
\(\lambda\)는 계수에 대해 어느 정도 벌점을 부여할 것인가를 결정한다. 이 값이 클수록 모델이 데이터에 오버피팅할 가능성이 낮아진다. lasso회귀 역시 이와 비슷한데, 다만 벌점 항에 유클리드 거리 대신 맨하탄 거리를 이용한다.
\[L_p = (\sum_{i = 1}^n|x_i|^p)^{\frac{1}{p}}\]
lambda [default = 1]
alpha [default = 0]
cf. Differences between L1 and L2 as Loss Function and Regularization
create_five_fold_cv <- function() {
set.seed(20190112)
lapply(createFolds(label, 5), function(idx) {
return(list(train = predictors[-idx, ],
validation = predictors[idx, ],
label_train = label[-idx],
label_val = label[idx]))
})
}
data_cv <- create_five_fold_cv()
params <- tibble(
eta = rep(c(.1, .5, .9), 3),
max_depth = rep(c(3, 6, 12), rep(3, 3))
)
xgb_avg_accuracy <- function(eta, max_depth) {
xgb_accuracy <- function(data) {
xgb <- xgboost(
data$train,
data$label_train,
params = list(
eta = eta,
max_depth = max_depth
),
objective = "binary:logistic",
nrounds = 100,
verbose = 0)
pred <- predict(xgb, data$validation)
mean(abs(data$label_val - pred) >= 0.5)
}
lapply(data_cv, xgb_accuracy)
}
avg_error <-
params %>%
mutate(
param_type = LETTERS[1:9],
avg_error = pmap(params, xgb_avg_accuracy)
) %>%
unnest(avg_error) %>%
mutate(avg_error = flatten_dbl(avg_error)) %>%
group_by(param_type) %>%
summarise(avg_error = mean(avg_error)) %>%
bind_cols(params, .) %>%
dplyr::select(-param_type) %>%
print()## # A tibble: 9 x 3
## eta max_depth avg_error
## <dbl> <dbl> <dbl>
## 1 0.1 3 0.329
## 2 0.5 3 0.334
## 3 0.9 3 0.343
## 4 0.1 6 0.331
## 5 0.5 6 0.355
## 6 0.9 6 0.375
## 7 0.1 12 0.347
## 8 0.5 12 0.374
## 9 0.9 12 0.385
XGBOOST HYPERPARAMETERS
xgboost 함수의 하이퍼파라미터들은 주로, 오버피팅, 정확도, 계산 복잡도 사이의 균형을 잡기 위해 사용된다.
비지도 학습은 여러 가지 서다 ㄷ른 목표들이 있을 수 있다. 어떤 경우에는 레이블이 정해진 응답변수가 없는는 상태에서 예측 규칙을 만드는 데 사용할 수 있다. 데이터의 의미 있는 그룹들을 찾기 위해 클러스터링을 사용할 수 있다. 예를 들면 웹사이트에서 사용자의 클릭 데이터와 인구통계 정보를 이용해 서로 다른 성격의 사요자들을 그룹화할 수 있을 것이다.
또 어떤 경우에는 데이터의 변수들을 관리할 수 있을 만한 수준으로 차원을 줄이는 것이 목표가 될 수도 있다. 이렇게 줄인 데이터는 회귀 혹은 분류 같은 예측 모델에 입력으로 사용할 수 있을 것이다.
마지막으로 변수와 레코드의 수가 아주 큰 상황이라면, 데이터 비지도 학습을 탐색적 데이터 분석의 연장으로 볼 수도 있다. 이때는 데이터와 다른 변수들 사이에 서로 어떤 관계가 있는지에 대한 통찰을 얻는 것이 목적이 된다. 비지도 학습은 변수들을 정밀하게 조사하고 관계를 밝혀내는 데 유용한 방법들을 제시한다.
Unsupervised Learning and Prediction
비지도 학습은 회귀와 분류에 대한 예측에 중요한 역할을 할 수 있다. 레이블이 없는 데이터에 대한 분류를 예측하기를 원할 수 있다. 예를 들어 위성 센서 데이터로부터 어떤 지역의 식생 유형을 예측하기를 원할 수 있다. 모델을 훈련시킬 반으변수가 없기 때문에 클러스터링은 지역을 분류하고 공통적인 패턴을 식별하기 위한 방법을 제공한다.
클러스터링은 특히 “cold-start problem”에 중요한 도구이다. 새로운 마케팅 홍보를 론칭하거나 잠재거인 새로운 형태의 사기나 스팸을 식별하는 캠페인에서는 모델을 훈련시킬 어떠한 반응변수도 가지고 있지 않다. 시간이 지나서 데이터가 모이면 시스템에 대해 좀 더 이해하고 전형적인 예측 모델을 학습할 수 있다. 하지만 클러스터링은 학습 프로세스를 좀 더 빨리 시작하도록 population segment를 식별함으로써 도와준다.
비지도 학습은 회귀와 분류 기술을 위한 block을 생성하는데에 중요하다. 빅 데이터를 가지고 만약 작은 subpopulation이 전반적인 population을 잘 대표하지 못한다면 훈련된 모델은 subpopulation에 대해 좋은 성능을 보이지 못한다. 클러스터링을 이용하여 subpopulation을 식별하고 식별하고 레이블을 지정할 수 있다. 그러면 서로 분리된 subpopualion들을 각각 다른 모델에 적합할 수 있다. 아니면 반대로, subpopulation 자체로 피쳐를 대표할 수 있어서 전체 모델이 subpopulation 정보를 명시적으로 예측변수로 고려하도록 하는 방법도 가능하다.
PCA은 수치형 변수가 어떻게 covery하는지 알아내는 방법이다. PCA는 선형판별분석의 비지도학습 버전이다.
oil_px <- sp500_data[, c("CVX", "XOM")]
pca <- princomp(oil_px)
pca$loadings##
## Loadings:
## Comp.1 Comp.2
## CVX 0.747 0.665
## XOM 0.665 -0.747
##
## Comp.1 Comp.2
## SS loadings 1.0 1.0
## Proportion Var 0.5 0.5
## Cumulative Var 0.5 1.0
첫 번째 주성분은 두 석유 회사 사이의 상관관계를 반영하는 CVX와 XOM의 평균을 의미한다. 반면 두 번째 성분은 CVX와 XOM의 주가가 달라지는 지점을 반영한다.
loadings <- pca$loadings
ggplot(oil_px, aes(CVX, XOM)) +
geom_point(alpha = .3) +
stat_ellipse(type = "norm", level = .99, size = 1) +
geom_abline(intercept = 0, slope = loadings[2, 1] / loadings[1, 1], linetype = "dashed", size = 1) +
geom_abline(intercept = 0, slope = loadings[2, 2] / loadings[1, 2], linetype = "dashed", size = 1) +
scale_x_continuous(breaks = -3:3, expand = c(0, 0)) +
scale_y_continuous(breaks = -3:3, expand = c(0, 0)) +
coord_fixed(xlim = c(-3, 3), ylim = c(-3, 3)) +
theme_bw() +
theme(
panel.grid.minor = element_line(color = "grey80", linetype = "dotted")
)In creating the first principal component, PCA arrives at the linear combination of predictor variables that maximizes the percent of total variance explained.
This linear combination then becomes the first “new” predictor, \(Z_1\)
PCA repeats this process, using the same variables, with different weights to create a second new predictor \(Z_2\). The weighting is done such that \(Z_1\) and \(Z_2\) are uncorrelated.
The process continues until you have as many new variables, or components, \(Z_i\) as original variables \(X_i\).
Choose to retain as many components as are needed to account for most of the variance.
The result so far is a set of weights for each compoent. The final step is to convert the original data into new principal component scores by applying the weights to the original values. These new scores and then be used as the reduced set of preditor variables.
syms <- c( 'AAPL', 'MSFT', 'CSCO', 'INTC', 'CVX', 'XOM', 'SLB', 'COP', 'JPM', 'WFC', 'USB', 'AXP', 'WMT', 'TGT', 'HD', 'COST')
top_sp <- sp500_data[sp500_data$X1 >= '2005-01-01', syms]
sp_pca <- princomp(top_sp)
screeplot(sp_pca)
Figure 6: A screeplot for a PCA of top stocks from the S&P 500
loadings <-
sp_pca$loadings[, 1:5] %>%
as.tibble() %>%
mutate(Symbol = rownames(sp_pca$loadings)) %>%
gather("Component", "Weight", -Symbol) %>%
mutate(positive = ifelse(Weight > 0, TRUE, FALSE))
# 책이랑 데이터가 다른 것 같다.
ggplot(loadings, aes(Symbol, Weight, fill = positive)) +
geom_bar(stat = "identity") +
facet_grid(Component ~ ., scales = "free_y") +
scale_fill_discrete(guide = FALSE) +
# guides(fill = FALSE) 특정 범례 제거
# theme(legend.position = "none") 모든 범례 제거
ylab("Component Loading") +
theme_bw() +
theme(axis.text.x = element_text(angle = -30, vjust = 1, hjust = 0))
Figure 7: The Loadings for the top five principal components of stock price returns
x와 y를 k = 4인 클러스터로 나눈다고 가정해보자.
\[\bar x_k = \frac{1}{n_k}\sum_{i \in Cluster k}x_i\]
\[\bar y_k = \frac{1}{n_k}\sum_{i \in Cluster k}y_i\] The sum of squares within a cluster is given by: \[SS_k = \sum_{i \in Cluster k}(x_i - \bar x_k)^2 + (y_i - \bar y_k)^2\] k-means는 \(SS_k\)를 최소가 되도록 레코드를 클러스터에 할당하는 방법이다.
df <- sp500_data[sp500_data$X1 >='2011-01-01', c('XOM', 'CVX')]
km <- kmeans(df, centers=4)
df$cluster <- factor(km$cluster)
centers <- data.frame(cluster = factor(1:4), km$centers)
ggplot(df, aes(XOM, CVX, color = cluster, shape = cluster)) +
geom_point(alpha = .3) +
geom_point(data = centers, aes(XOM, CVX), size = 3, stroke = 3) +
theme_bw()
Figure 8: The clusters of K-means applied to stock price data for ExxonMobil and Chevron(the two cluster centers in the dense area are hard to distinguish
정확한 k-means해는 계산하기 어려우므로, heuristic algorithm
사용자가 미리 정해준 K값과 클러스터 평균의 초깃값을 가지고 알고리즘을 시작하며, 아래 과정을 반복한다.
알고리즘은 각 레코드에 대한 클러스터 할당이 더 이상 변화가 없을 때 수렴한다. 첫 번째 단계에서 클러스터 평균의 초깃값을 설정할 필요가 있다. 보통은 각 레코드를 K개의 클러스터들 가운데 하나에 랜덤하게 할당한 후 그렇게 얻은 클러스터들의 평균을 사용한다.
이 방법이 항상 최적의 답을 준다는 보장이 없기 때문에, 랜덤하게 초깃값을 변화시켜가며 알고리즘을 여러 번 돌려봐야 한다. 일단 초깃값이 정해지게 되면 반복 루프를 통해 K-menas는 클러스터 내 제곱합이 최소가 되도록 하는 해를 얻게 된다.
R에서 kmeans함수의 nstart 인자를 이용해 랜덤하게 초깃값을 다르게 설정해 알고리즘을 시행할 횟수를 설정할 수 있다. 다음은 5개의 클러스터를 찾기 위해 K평균 알고리즘을 서로 다른 초깃값을 이용해 10번 수행하는 코드이다.
syms <- c( 'AAPL', 'MSFT', 'CSCO', 'INTC', 'CVX', 'XOM', 'SLB', 'COP',
'JPM', 'WFC', 'USB', 'AXP', 'WMT', 'TGT', 'HD', 'COST')
df <- sp500_data[sp500_data$X1 >= '2011-01-01', syms]
km <- kmeans(df, centers=5, nstart=10)10번의 시도 가운데 가장 좋은 성능을 보인 결과를 볼 수 있다. iter.max변수를 이용해 각 초기 설정별 알고리즘의 최대 반복 횟수를 설정할 수 있다.
kmeans함수에서 가장 중요한 두 아웃풋은 클러스터의 크기와 클러스터 평균이다.
km$size## [1] 285 106 288 266 186
km$centers %>%
as.tibble() %>%
gather("Symbol", "Mean") %>%
mutate(cluster = paste("Cluster", rep(1:5, length(as.data.frame(km$centers)))),
is_positive = Mean > 0) %>%
ggplot(aes(Symbol, Mean, fill = is_positive)) +
geom_bar(stat = "identity") +
facet_grid(cluster ~., scale = "free_y") +
ylab("Component Loading") +
theme_bw()
Figure 9: The means of the variables in each cluster(cluster means)
K 평균 알고리즘을 사용하려면 클러스터의 개수를 지정해야한다. 클러스터 개수는 응용 분야에 따라 결정되기도 한다. 고객들을 분류할 때 관리를 위 한 고려 사항에 따라 고객들을 분류할 그룹의 개수를 미리 결정할 수 있다.
실무 혹은 관리상의 고사 사항에 따라 클러스터 개수를 미리 결정하기가 어려울 경우, 통계적 접근 방식을 사용할 수 있다.
elbow method은 언제 클러스터 세트가 데이터의 분산의 ’대부분’을 설명하는지를 알려준다. CV를 사용하여 부분적으로 평가할 수 있다.
kmean_betweenss <- function(k) {
cluster <- kmeans(df, k, nstart = 50, iter.max = 100)
cluster$betweenss
}
# total all SS
tot_all_SS <- kmeans(df, 14, nstart = 50, iter.max = 100)$totss
# Set maximum cluster
max_k <- 14
# Run algorithm over a range of k
bss <- sapply(2:max_k, kmean_betweenss)
# Create a data frame to plot the graph
elbow <- tibble(num_clusters = 2:max_k, var_explained = bss / tot_all_SS)
# Plot the graph
ggplot(elbow, aes(num_clusters, var_explained)) +
geom_point() +
geom_line() +
scale_x_continuous(breaks = seq(1, max_k, by = 1)) +
labs(x = "% Variance Explained", y = "Number of Clusters") +
theme_bw()
Figure 10: The elbow method applied to the stock data
클러스터 개수를 결정하는 데에는 통계 이론이나 정보이론에 바탕을 둔 좀 더 형식적인 방법들도 있다. 예를 들어 “gap” statistic을 사용하여 클러스터 개수를 결정하는 데에 이용하기도 한다. 하지만 이러한 이론적 접근법은 대부분의 응용 분야에서 필요하지 않고, 심지어 적합하지 않을 수도 있다.
The 5 Clustering Algorithms Data Scientists Need to Know
K-means보다 유연하고, 수치형 변수가 아니어도 쉽게 적용이 가능하다. 이상치나 비정상적이 그룹 또는 레코드를 발견하는 데에 좀 더 민감하다. 계층적 클러스터링은 또한 직관적인 시각화가 가능하여 클러스터를 해석하기가 수월하다.
계층적 클러스터링의 유연성에는 비용이 따른다. 다시 말해서 계층적 클러스터링은 수백만 개의 레코드가 있는 데이터에는 적용할 수 없다. 수만 개의 레코드로 이루어진 적당한 크기의 데이터의 경우에도 계층적 클러스터링을 위해서는 상대적으로 많은 computing resource가 필요할 수 있다. 실제로, 계층적 클러스터링은 대부분 상대적으로 데이터 크기가 작은 문제에 적용된다.
아래 두 가지 기본 구성 요소를 기반으로 한다.
데이터 자체보다 쌍 거리 \(d_{i,j}\)에 따라 동작한다.
syms1 <- c('GOOGL', 'AMZN', 'AAPL', 'MSFT', 'CSCO', 'INTC', 'CVX',
'XOM', 'SLB', 'COP', 'JPM', 'WFC', 'USB', 'AXP', 'WMT', 'TGT', 'HD', 'COST')
# take transpose: to cluster companies, we need the stocks along the rows
df <- t(sp500_data[sp500_data$X1 >= '2011-01-01', syms1])
d <- dist(df)
hcl <- hclust(d)plot(hcl)나무의 잎은 레코드에 대응한다. 나뭇가지의 길이는 클러스터 간의 비유사성 정도를 나타낸다. 구글과 아마존에 대한 수익률은 다른 주식대 대한 수익률과 상당히 다른 것을 볼 수 있다. 다른 주식들을 자연스럽게 그룹으로 형성된다. 에너지, 금융, 소매기업들이 각각 하위 트리로 구분된 것을 볼 수 있다.
# 클러스터의 구체적인 수
cutree(hcl, k = 4)## GOOGL AMZN AAPL MSFT CSCO INTC CVX XOM SLB COP JPM WFC
## 1 2 3 3 3 3 4 4 4 4 3 3
## USB AXP WMT TGT HD COST
## 3 3 3 3 3 3
구글과 아마존은 각자 자신만의 독자적인 클러스터에 속해 있음을 볼 수 있다. 석유 주식((XOM, CVS, SLB, COP))은 모두 한 클러스터에 속한다. 마지막으로 나머지 주식들이 남은 하나의 클러스터에 속한다.
거리로 유클리드 거리를 사용한다고 하자. 비유사도를 측정하는 한 가지 방법은 A와 B 사이의 모든 레코드 쌍의 최대 거리를 사용하는 complete linkage
\[D(A, B) = max\space d(a_i, b_j)\space for\space all\space pairs\space i, j\] 1. 데이터의 모든 레코드에 대해, 단일 레코드로만 구성된 클러스터들로 초기 클러스터 집합을 만든다. 2. 모든 쌍의 클러스터 \(k, l\)사이의 비유사도 \(D(C_k, C_l)\)을 계산한다. 3. \(D(C_k, C_l)\)에 따라 가장 가까운 두 클러스터 \(C_k\)와 \(C_l\)을 병합한다. 4. 둘 이상의 클러스터가 남아 있으면 2단계로 다돌 돌아간다. 그렇지 않고 클러스터가 하나 남는다면 알고리즘을 멈춘다.
비유사도를 측정하는 네 가지 일반적인 지표는 complete linkage, single linkage, average linkage, and minimum variance single linkage방법은 두 클러스터에 레코드 사이의 거리를 최소화 한다. \[D(A, B) = min\space d(a_i, b_j)\space for\space all\space pairs\space i, j\] 이는 greedy방법이며, 결과로 나온 클러스터는 서로 크게 다른 요소들을 포함하는 일도 생길 수 있다. average linkage method는 모든 거리 쌍의 평균을 사용하는 방법으로 이는 single linkage method와 complete linkage method사이를 절충한 방법이다. 마지막으로 Ward’s method이라고도 하는 최소분산 방법은 클러스터 내의 제곱합을 최소화하므로 k-means와 유사하다고 할 수 있다.
df <- sp500_data[sp500_data$X1 >='2011-01-01', c('XOM', 'CVX')]
d <- dist(df)
hclust_single <- cutree(hclust(d, method = "single"), k = 4)
hclust_avg <- cutree(hclust(d, method = "average"), k = 4)
hclust_complete <- cutree(hclust(d, method = "complete"), k = 4)
hclust_ward <- cutree(hclust(d, method = "ward.D"), k = 4)
df <-
df %>%
mutate(
hclust_single = hclust_single,
hclust_avg = hclust_avg,
hclust_complete = hclust_complete,
hclust_ward = hclust_ward
) %>%
gather("cluster_method", "cluster_class", -XOM, -CVX) %>%
mutate(cluster_method = factor(cluster_method, labels = c("average", "complete", "single", "ward.D")),
cluster_class = as.factor(cluster_class))
ggplot(df, aes(XOM, CVX, color = cluster_class, shape = cluster_class)) +
geom_point(alpha = .3) +
guides(color = guide_legend(override.aes = list(alpha = 1),
title = "cluster"),
shape = guide_legend(title = "cluster")) +
# labs(color = "cluster", shape = "cluster") + # color shape 둘 다 바꿔야한다!!
facet_wrap(~cluster_method) +
theme_bw()
Figure 11: A comparison of measures of dissimilarity applied to stock data
single linkage는 거의 모든 점을 하나의 클러스터에 할당했다. 최소분산 방법(Ward.D)을 제외한 다른 방법은 외곽에 요소가 몇 없는 특이값들로 이루어진 작은 클러스터를 만들었다. 최소분산 방법은 k-means와 가장 유사한 결과를 보인다.
예를 들어, 전반적으로는 서로 비슷하지만 모든 데이터가 반드시 서로 가까울 필요는 없는 그룹(예를 들어 수익의 분산이 큰 기술 주식들)과 서로 비슷하면서 데이터들이 아주 가까이에 있는 또 다른 그룹(예를 들어 분산이 적은 유틸리티 주식)이 함께 있는 경우에 사용할 수 있다.
모델 기반 크러스터링의 핵심 아이디어는 각 레코드가 K개의 다변량정규분포 중 하나로부터 발생했다고 가정하는 것이다. 여기서 K가 클러스터의 개수를 의미한다. 각 분포는 서로 다른 평균 \(\mu\)와 공분산행렬 \(\Sigma\)를 갖는다. 예를 들어 X와 Y라는 두 개의 변수가 있는 경우 각 행\(X_i, Y_i\)은 K개의 분포 \(N_2(\mu_1, \Sigma_1),N_2(\mu_2, \Sigma_2), ..., N_2(\mu_K, \Sigma_K)\) 중 하나에서 샘플링된 것으로 모델링된다.
library(mclust)## Warning: package 'mclust' was built under R version 3.5.2
## Package 'mclust' version 5.4.2
## Type 'citation("mclust")' for citing this R package in publications.
##
## Attaching package: 'mclust'
## The following object is masked from 'package:purrr':
##
## map
df <- sp500_data[sp500_data$X1>='2011-01-01', c('XOM', 'CVX')]
mcl <- Mclust(df)
summary(mcl)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VEE (ellipsoidal, equal shape and orientation) model with 2
## components:
##
## log.likelihood n df BIC ICL
## -2255.125 1131 9 -4573.528 -5075.657
##
## Clustering table:
## 1 2
## 168 963
cluster <- as.factor(predict(mcl)$classification)
ggplot(df, aes(XOM, CVX, color = cluster)) +
geom_point(size = 0.75) +
theme_bw()
Figure 12: Two clusters are obtained for stock return data using mclust
summary(mcl, parameters = TRUE)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VEE (ellipsoidal, equal shape and orientation) model with 2
## components:
##
## log.likelihood n df BIC ICL
## -2255.125 1131 9 -4573.528 -5075.657
##
## Clustering table:
## 1 2
## 168 963
##
## Mixing probabilities:
## 1 2
## 0.2739885 0.7260115
##
## Means:
## [,1] [,2]
## XOM -0.04362218 0.05792282
## CVX -0.21109525 0.07375447
##
## Variances:
## [,,1]
## XOM CVX
## XOM 1.044671 1.065190
## CVX 1.065190 1.912748
## [,,2]
## XOM CVX
## XOM 0.2998935 0.3057838
## CVX 0.3057838 0.5490920
predict.Mclust, summary.Mclust 참고
플랏을 보면 주식 데이터는 정규분포 모양을 갖고 있는 것 처럼 보인다. 사실 더 정확히는, 주식 수익률은 정규분포보다 긴 꼬리 분포를 따른다고 볼 수 있다. 이를 처리하기 위해 mclust는 데이터의 대부분에 대한 분포를 피팅하려고 하고 그 결과 두 번째 분포가 더 큰 분산을 갖도록 피팅된다.
k-means와 hierachical 클러스터링과 달리 mclust는 클러스터 수(이 경우 2)를 자동으로 선택한다. 바로 Bayesian information criteria(BIC)값이 가장 큰 클러스터의 개수를 선택하도록 동작하기 때문이다.
plot(mcl, what = "BIC", ask = FALSE)
Figure 13: BIC scores for the stock return data for different numbers of clusters(components)
k-means의 클러스터 개수를 선택하는데 사용한 elbow plot과 유사하지만, 여기서는 설명되는 분산 비율 대신 BIC가 사용되었다. 한 가지 큰 차이점은 선이 하나가 아니라, mclust 함수가 무려 14개의 다른 선을 동시에 보여준다는 것이다! 이는 mclust가 실제로 각 클러스터 크기에 대해 14개의 다른 모델을 피팅했고 궁극적으로 가장 적합한 모델을 선택했다는 것을 잘 보여준다.
모델 기반 클러스터링에는 몇 가지 한계가 있다. 이 방법은 기본적으로 데이터들이 모델을 따른다는 가정이 필요하며, 클러스터링 결과는 이 가정에 따라 매우 다르다. 필요한 계산량 역시 계층적 클러스터링보다 높으므로 대용량 데이터로 확장하기가 어렵다. 마지막으로 알고리즘이 다른 방법들보다 더 복잡하고 이용하기가 어렵다.
데이터의 크기가 조정되지 않으면 PCA, K-means, 혹은 기타 클러스터링 방법은 큰 값을 갖는 변수들에 의해 결과가 좌우되고 작은 값을 갖는 변수들은 무시된다.
범주형 데이터는 일부 클러스터링 과정에서 특별한 문제를 일으킬 수 있다. KNN에서와 마찬가지로, 순서가 없는 요인변수는 일반적으로 one-hot encoding을 사용하여 이진(0/1) 변수 집합으로 변 환한다. 이러한 이진변수는 다른 데이터와 스케일이 다를 뿐만 아니라 PCA와 K-means 같은 기법을 사용할 때 이진변수가 두 가지의 값만 가질 수 있다는 것 때문에 문제가 될 수 있다.
PCA에서 상관행렬을 사용하면 z점수를 사용하는 것과 같은 효과이다.
변수들이 서로 동일한 규모로 측정되고 상대적 중요성을 정확하게 반영하는 경우(예를 들어 주가 변동)조차도 변수의 스케일을 재조정하는 것이 유용할 수 있다.
syms <- c('AMZN', 'GOOGL', 'AAPL', 'MSFT', 'CSCO', 'INTC', 'CVX', 'XOM', 'SLB', 'COP', 'JPM', 'WFC', 'USB', 'AXP', 'WMT', 'TGT', 'HD', 'COST')
top_sp1 <- sp500_data[sp500_data$X1 >='2005-01-01', syms]
sp_pca1 <- princomp(top_sp1)
screeplot(sp_pca1)
Figure 14: A screeplot for a PCA of top stocks from the S&P 500 including GOOGL and AMZN
round(sp_pca1$loadings[, 1:2], 3)## Comp.1 Comp.2
## AMZN 0.445 0.874
## GOOGL 0.857 -0.478
## AAPL 0.072 0.021
## MSFT 0.036 0.006
## CSCO 0.029 0.003
## INTC 0.027 0.006
## CVX 0.090 0.037
## XOM 0.080 0.021
## SLB 0.110 0.030
## COP 0.058 0.024
## JPM 0.071 0.009
## WFC 0.053 0.009
## USB 0.042 0.006
## AXP 0.079 0.024
## WMT 0.040 0.007
## TGT 0.064 0.025
## HD 0.051 0.033
## COST 0.071 0.034
처음 두 가지 주성분이 거의 완전히 GOOGL과 AMZN에 의해 지배되고 있다. 이는 GOOGL과 AMZN의 주가 움직임이 전체 변동성의 대부분을 지배하기 때문이다.
이러한 상황에서는 변수를 스케일링 해서 포함하거나 dominant 변수를 분석으로부터 제외하고 별도로 처리할 수 있다. 어떤 방법이 항상 옳다고는 볼 수 없으며 응용 분야에 따라 달라진다.
범주형 데이터가 있는 경우에는 순서형 변수 또는 이진형(더미) 변수를 사용하여 수치형 데이터로 변환해야 한다. 데이터를 구성하는 변수들에 연속형과 이진형 변수가 섞여 있는 경우에는 비슷한 스케일이 되도록 변수의 크기를 조정해야 한다. 이를 위한 대표적인 방법은 Gower’s Distance를 사용하것 ㄱ이다.
고워 거리의 기본 아이디어는 각 변수의 데이터 유형에 따라 거리 지표를 다르게 적용하는 것이다.
고워 거리는 다음과 같이 계산한다.
x <- defaults[1:5, c('dti', 'payment_inc_ratio', 'home', 'purpose')]
x# A tibble: 5 × 4
dti payment_inc_ratio home purpose
<dbl> <dbl> <fctr> <fctr>
1 1.00 2.39320 RENT car
2 5.55 4.57170 OWN small_business
3 18.08 9.71600 RENT other
4 10.08 12.21520 RENT debt_consolidation
5 7.06 3.90888 RENT other
cluster패키지의 daisy함수를 사용하면 고워 거리를 계산할 수 있다.
library(cluster)
daisy(x, metric='gower')Dissimilarities :
1 2 3 4
2 0.6220479
3 0.6863877 0.8143398
4 0.6329040 0.7608561 0.4307083
5 0.3772789 0.5389727 0.3091088 0.5056250
거리가 모두 0과 1사이인 것을 볼 수 있다. 가장 거리가 먼 레코드 쌍은 2번과 3번 레코드이다. 이 둘은 home_변수나 purpose 변수에 대해서도 값이 다르며 dti 변수와 payment_inc_ratio 변수 역시 차이가 크다. 반면 home_변수나 purpose 변수가 동일한 값을 갖는 레코드 3번과 5번은 거리가 가장 적다.
daisy함수의 거리 행렬 결과에 hclust 함수를 이용하여 계층적 클러스터링을 적용해볼 수 있다.
df <- defaults[sample(nrow(defaults), 250),
c('dti', 'payment_inc_ratio', 'home', 'purpose')]
d = daisy(df, metric='gower')
hcl <- hclust(d)
dnd <- as.dendrogram(hcl)
plot(dnd, leaflab='none')다음 코드를 사용하면 하위 트리 중 하나(cut을 사용해 0.5를 기준으로 잘랐을 때 가장 왼쪽에 있는 그룹)에 포함된 레코드들을 살펴볼 있다.
dnd_cut <- cut(dnd, h = .5)
df[labels(dnd_cut$lower[[1]]),]# A tibble: 9 × 4
dti payment_inc_ratio home purpose
<dbl> <dbl> <fctr> <fctr>
1 24.57 0.83550 RENT other
2 34.95 5.02763 RENT other
3 1.51 2.97784 RENT other
4 8.73 14.42070 RENT other
5 12.05 9.96750 RENT other
6 10.15 11.43180 RENT other
7 19.61 14.04420 RENT other
8 20.92 6.90123 RENT other
9 22.49 9.36000 RENT other
이 하위 트리에 속한 레코드들은 모두 대출 목적이 other이고 주거 형태가 RENT이다. 모든 하위 트리에서 이러한 형태의 구분이 이루어지는 것은 아니지만 이는 범주형 변수 값이 비슷한 데이터들이 한 클러스터로 그룹화되는 경향이 있다는 것을 보여준다.
k-means와 PCA는 연속형 변수에 가장 적합하다. 데이터 집합의 크기가 더 작아질수록 고워 거리를 사용하여 계층적 클러스터링을 하는 것이 좋다. 원친적으로는 이진형 혹은 범주형 데이터에도 K-means를 적용할 수 있다. 범주형 데이터의 경우에는 일반적으로 one-hot encoding 방법을 이용해 수치형으로 변환할 수 있다. 하지만 실무에서는 k-means와 PCA를 이진형 데이터와 함께 사용하는 것은 어려울 있다.
표준 z점수를 사용할 경우, 이진형 변수가 클러스터 결과에 지대한 영향을 미치게 된다. 0/1변수의 경우, 0 또는 1의 값인 레코드를 모두 한 클러스터에 포함되므로 K-means평균에서 클러스터내 제곱합이 작아지기 때문이다.
df <- model.matrix(~ -1 + dti + payment_inc_ratio + home + pub_rec_zero, data=defaults)
df0 <- scale(df)
km0 <- kmeans(df0, centers=4, nstart=10)
centers0 <-scale(km0$centers, center=FALSE,
scale=1/attr(df0, 'scaled:scale'))
scale(centers0, center=-attr(df0, 'scaled:center'), scale=F)dti payment_inc_ratio homeMORTGAGE homeOWN homeRENT pub_rec_zero
1 17.02 9.10 0.00 0 1.00 1.00
2 17.47 8.43 1.00 0 0.00 1.00
3 17.23 9.28 0.00 1 0.00 0.92
4 16.50 8.09 0.52 0 0.48 0.00
상위 4개 클러스터가 요인변수들과 밀접한 관련이 있다는 것을 볼 수 있다. 이러한 문제를 피하기 위해 이진형 변수의 크기를 다른 변수들보다 작은 값으로 조정할 수 있다. 이러한 문제를 피하기 위해 이진형 변수의 크기를 다른 변수들보다 작은 값으로 조정할 수 있다. 아니면 데이터의 그기가 아주 큰 경우에는 특정 범주 값들에 따라 서로 다른 하위 집합에 클러스터링을 적용할 수도 있다. 예를 들어, 주택 담보 대출이 있는 사람인지, 주택을 완전히 소유한 사람인지, 주택을 임대하고 있는 사람인지에 따라 대출 데이터를 쪼개서 각각의 하위 그룹에 클러스터링을 개별적으로 적용할 수 있다.
주성분분석와 k-means 클러스터링은 수치형 데이터의 차원을 축소하기 위해 주로 사용되는 방법들이다. 의미 있는 데이터 축소를 보장하기 위해서는 데이터의 스케일을 적절히 조정해야 한다.
그룹들 간의 구분이 분명하고 고도로 구조화된 데이터를 클러스터링할 경우, 어떤 방법을 사용하든 결과는 비슷할 가능성이 높다. 물론 방법마다 장점이 있다. K-means은 매우 큰 데이터로 확장이 가능하고 이해하기 쉽다. 계층적 클러스터링은 수치형과 범주형이 혼합된 데이터 유형에 적용이 가능하며 직관적인 시각화 방법이 존재한다. 모델 기반 클러스터링은 휴리스틱한 방법들과 달리 통계 이론에 기초를 두고 있으며 더 엄밀한 접근 방식을 제시한다. 그러나 데이터가 커지면 K-means가 가장 많이 사용되는 방법이다.
대출이나 주식 데이터, 그리고 데이터 과학자가 직면할 대다수 데이터는 노이즈가 많다. 이런 경우, 사용 기법에 따라 결과에 극명한 차이를 가져온다. K-means, 계층적 클러스터링, 그리고 특히 모델 기반 클러스터링은 모두 매우 다른 솔루션을 생성한다. 데이터 과학자는 이럴 때 어떻게 해야 할까? 불행하게도 선택을 돕는 간단한 법칙 따위는 없다. 궁극적으로 데이터 크기나 응용 분야의 목표에 따라 사용되는 방법은 달라지게 된다.